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These lecture notes cover the statics and glassy dynamics of granular media. Most of the lectures were in fact devoted to 
'force propagation' models. We discuss the experimental and theoretical motivations for these approaches, and their conceptual 
connections with Edwards' thermodynamical analogy. One of the distinctive feature of granular media (common to many other 
'jammed' systems) is indeed the large number of metastable states that are macroscopically equivalent. We present in detail 
the (scalar) g-model and its tensorial generalization, that aim at modelling the existence of force chains and arching effects 
without introducing any displacement field. The contrast between the hyperbolic equations obtained within this line of thought 
and elliptic (elastic) equations is emphasized. The role of disorder on these hyperbolic equations is studied in details using 
perturbative and diagrammatic methods. Recent (strong disorder) force chain network models are reviewed, and compared 
with the experimental determination of the force 'response function' in granular materials. We briefly discuss several issues 
(such as isostaticity and generic marginality) and open problems. At the end of these notes, we also discuss the basic dynamical 
properties of weakly tapped granular assemblies, and stress the phenomenological analogies with other glassy materials. Simple 
models that account for slow compaction and dynamical heterogeneities are presented, that are inspired by 'free-volume' ideas 
and Edwards' assumption. A connection with the theory of fluctuating random surfaces, also noted recently by Castillo et al., 
is suggested. Finally, we discuss how the 'trap model' can be adapted to granular materials, such that more subtle 'memory' 
effects can be accounted for. 



I. INTRODUCTION 
A. Basic phenomenology 

Although granular materials are made of classical particles of macroscopic size, they exhibit a host of interesting and 
sometimes counter-intuitive properties, which are of interest both to the academic community and for industrial applications: 
enormous amounts of 'granular assemblies' are routinely handled (stored, transported, mixed together, etc.) The reason why 
these systems are still not fully understood yet is that they require a proper statistical treatment of collective effects. Although 
the physics at the grain level is reasonably well understood, the behaviour of a large assembly of these grains, with strongly 
non-linear interactions, demand concepts and methods that only recently emerged from the study of disordered systems in 
statistical mechanics. The main property of granular materials, that leads both to most of the non trivial phenomenology 
and to most of the theoretical difficulties, is the existence of a large number of different microscopic metastable states that are 
macroscopically equivalent. It turns out, quite interestingly, that this feature (metastability) is common to a wider class of 
materials in their 'jammed' state. This includes glasses, colloids, compressed emulsions and foams, spin-glasses, vortex glasses 
and other collectively pinned structures, etc. 

There is a vast body of experimental results on granular materials that we do not aim to cover here (for reviews, see 
[81,80,56,57]). We will mainly restrict to dry grains in static and weakly driven situations. Correspondingly, we will focus 
on the statistical properties of the static states (and in particular the distribution of stresses) and on the glassy dynamics of 
gently 'tapped' assemblies, that slowly evolve from one static state to another. We have therefore excluded from these lectures, 
because of lack of time, 'strongly' driven situations, such as granular flows, avalanches and surface flows, dune formation, etc. 
This is not to suggest that these problems are less interesting and that collective effects are irrelevant. Quite on the contrary, 
the dynamics of strongly driven granular systems also displays remarkable effects, such as coUisional clustering which generates 
non trivial spatial structures in granular flows, and invalidates simple hydrodynamical descriptions [72,96,16]. Some recent 
papers on these matters can be found in references [80,57,112,98,113]. 

Stress patterns in dry granular media exhibit some rather unusual features when compared to either liquids or elastic solids. 
For example, the vertical pressure below conical sand-piles does not follow the height of material above a particular point. 
Depending on the way the pile is prepared, it shows a minimum underneath the apex of the pile [128,29,140] when the pile is 
built from a point source, and a broad, flat maximum when it is built layer by layer from a uniform 'rain' of grains. Furthermore, 
local stress fluctuations are large, sometimes on length scales much larger than the grain size. For example, repeatedly pouring 
the very same amount of powder in a silo results in fluctuations of the weight supported by the bottom plate of 20% or more 
[30,141]. Weak perturbations of the packing can sometimes cause large rearrangements [44,50]. Qualitatively, these features 
are attributed to the presence of stress paths which can focus the stress fleld into localized regions and also deflect it to cause 
"arching" (see [49] for early qualitative experiments). More quantitative experiments were performed in [93,29,13], where the 
local fluctuations of the normal stress deep inside a silo or at the base of a sandpile were measured. It was found that the stress 
probability distribution is rather broad, decaying exponentially for large stresses. This behaviour was also found in numerical 
simulations [114], and more recently in other situations [109], such as compressed emulsions [31]. Similarly, standard 'triaxial' 
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test experiments used to determine the elastic properties of materials from the (macroscopic) relation between stresses and 
deformations show highly irreproducible, hysteretic behaviour which only seem to converge towards a well defined curve after 
a large number of deformation cycles have been imposed to the granular system in order to 'anneal' it down to a reproducible 
state. 

The dynamics of slowly driven granular systems also exhibit unusual features when compared to either liquids or solids, 
and has actually much in common with glasses. In particular, the way these systems very slowly compact when vibrated, the 
unusual dependence of the density on the system history, etc. has strong similarities with the properties of glassy materials. At 
the phenomenological level, the dynamics of these systems appears to be a succession of hops between diflferent (metastable) 
equilibrium states. The understanding of static and weakly driven granular assemblies therefore require a proper description 
of the statistical properties of these 'blocked' configurations. We now discuss these issues in the perspective of the present 
lectures. ^ 



B. Theoretical issues 

1. Static properties 

How can one then describe the statics of granular materials on large length scales? The basic problem stems from the fact 
that the equilibrium equations for the stress tensor are not sufficient to determine the stress. For example, in two dimensions the 
stress tensor has three independent components, but there arc only two equilibrium equations. Some additional assumptions 
about the properties of the material must be provided. For example, the assumption that the material is elastic and follows 
Hooke's law gives extra constraints on the stress tensor and allows one to solve the static problem as soon as some appropriate 
boundary conditions are given. For granular materials, the standard procedure is to use elastic or elasto-plastic theories from 
soil mechanics [148]. However, the relation between force chains on short length scales and an elasto-plastic description on 
large length scales is far from obvious. To our knowledge, no systematic procedure has ever been proposed to go from the 
mechanics at the grain level to coarse-grained equations that would justify the use of an elasto-plastic framework and estimate 
the parameters of the theory (effective elastic moduli, etc.). One of the main difficulty is that some indeterminacy exists already 
at the grain level, since many different configurations of the contact forces are allowed and satisfy local equilibrium. This leads 
to several conceptual problems: even if an elastic-like description of small perturbations around an arbitrary reference state 
(such as sound waves, for example) might make sense in general, the description of - say a conical sandpile using elasto- 
plasticity theory [33] requires the identification of a (zero stress) reference state from which deformations can be defined, at 
least for some regions of the pile. In the case of a pile of infinitely hard grains (which should be the correct benchmark for an 
assembly of grains with a Young modulus much larger than the gravity induced stresses) that rests in one particular metastable 
state (among a large rmmber of macroscopically equivalent ones), switching the gravity back to zero will hardly affect the 
packing. Each of these metastable states can thus equally well be taken as a reference state; on the other hand, it is precisely 
the stress pattern in one of these 'native' metastable state (i.e. obtained when grains come to rest without further tapping) 
that one wants to predict. One peculiarity of (dry) granular materials is the absence of tensile stresses between the grains; the 
cohesion of the assembly is therefore induced by the applied stress itself and the zero stress state is ill defined. 

Even the description of small perturbations around a given reference state might be problematic. For example, the existence 
of a (large volume) limiting curve relating incremental stresses and deformations requires, as already mentioned above, at 
least some 'annealing' procedure to define a reproducible initial state. This limiting curve might not even exist in the absence 
of firiction [44]. Even for moderate deformations, following the so-called consolidation phase, the response to cyclic loads in 
standard triaxial tests shows some significant irreversibility. 

The absence of any obvious deformation field from which the stress tensor may be constructed has motivated an alternative, 
'stress-only' approach [18,145,21,42]. The basic tenet of these theories is that in equilibrium, some (history dependent) large 
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scale relations between the components of the stress tensor should exist. These relations should be determined by the global 
statistical features of the particular metastable state in which the packing sits but not on its microscopic details, nor on the 
particular loading conditions, provided these do not lead to further rearrangements of the packing. Much as random collisions 
between molecules give rise, on large length scale, to well defined hydrodynamical equations, the hope is that an appropriate 
coarse-graining of the local force balance equations leads, on large length scales, to the missing 'closure' equation that allows 
to solve for the static equilibrium. (For rather formal attempts in this direction, see [62,2].) A well known relation of this 
type arises from the assumption that the material is everywhere on the verge of plastic failure, leading to a Mohr-Coulomb 
(non- linear) relation between the stress components [106], but we will motivate and discuss simpler relations below, based both 
on symmetry arguments [18] and on the consideration of simple rules for the transfer of stresses between adjacent grains [41]. 
The consequence of a fixed relation between the components of the stress tensor is that stresses obey an hyperbolic equation, as 
compared to the elliptic equations encountered in elasticity theory. This means that stresses 'propagate' or arc 'transmitted' 
along lines: as discussed below, the characteristics of this hyperbolic equation are the mathematical transcription of the force 
chains that are well known to exist in granular materials [35] . 



2. Tapping and non thermal ensembles 

As mentioned above, many different packings and configurations of the contact forces are compatible with the local equilibrium 
of each grain for a given macroscopic situation. This is actually intimately related to the fact that stresses in granular media 
often show large fluctuations; some kind of averaging is therefore needed to obtain reproducible results. In the case of sand- 
piles, one must repeat the construction of the pile several times, and use a pressure gauge that averages over a sufficiently 
large number of grains, in order to obtain a satisfactory stress profile that a statistical theory of blocked states should predict. 
Another possibility is to vibrate the packing such as to make it probe, during its evolution, several equilibrium states with the 
same macroscopic geometry. The natural question is then: with which statistical weight the different equilibrium (blocked) 
states appear in a given experiment? To what extent are these weights dependent on the dynamics that leads to the blocked 
states? Is the ensemble of 'native' (as-built) packings identical to the ensemble of packings obtained under tapping? 

The simplest answer, proposed more than ten years ago by Sam Edwards, is to postulate that all blocked states with a given 
density are equiprobable [60]. This micro-canonical assumption defines what is now called the 'Edwards ensemble' [86]; it turns 
out that several toy models of jammed systems do obey, either exactly or to a good approximation, Edwards' prescription 
[5,52,127,45]. This is a first step towards a 'thermodynamical' description of out-of-equilibrium, dissipative systems [3,97,8,4]. 
However, several remarks of various nature should be made here. 

• First, the analogy between tapping strength and temperature. In many cases, this is a useful intuitive guide and several 

experiments discussed in section X do indeed confirm the phenomcnological analogies between the two. However, tapping 
is a long-wavelength excitation, whereas temperature in solid state physics is thought to give rise to very short wavelength 
fluctuations. Although the long-wavelength excitation probably cascades down, through collisions, to short wavelengths, 
ideas such as detailed balance and activated processes might be affected by the correlated nature of the noise. In this 
respect, the non trivial clustering patterns induced by the dissipative collisions might also obliterate simple ideas on the 
statistics of blocked states. 

• One must distinguish at least two types of tapping excitations. One would be very gentle taps, that are insufficient to 

change the packing geometry, but do change the contact forces for each grain. In this case, tapping induces a random 
walk in 'force space', but for a fixed configuration of the grains. The Edwards hypothesis in this restricted case is to 
assign a uniform weight for all force configurations that (i) lead to static equilibrium (forces and torques on each grain 
add up to zero) and (ii) satisfy the Coulomb inequality at each contact. One can also drive the system with an amplitude 
such that the motion of grains is possible, in which case the dynamics is a random walk both in force space and in packing 
space. The extended Edwards ensemble in this case is to assign equal weight to any packing and any force configuration 
such that equilibrium is obeyed and the Coulomb inequality satisfied. In principle, the Edwards prescription could be 
correct in one case and not in the other, or in both, or in none, or more complicated situations still. 

• The Edwards prescription is however ambiguous for continuous variables. In the 'gentle' tap case, one is tempted to 
interpret Edwards' measure as follows. Let us call /" the contact force on the Q-th contact of the j-th grain, and 
the position of the contact point. We call ii the friction coefficient, and the indices A'^ and T refer to the normal and 
tangential components of the force. The natural Edwards measure reads: 



J2fns(j2f^x^^]l[®H"N-\fr,T\) 



(1) 



which is a formal way to impose the constraints on each grain. (© is the step function). However, this assumes that the 
a priori measure on the forces is uniform, which is reasonable but not obvious. The usual microcanonical ensemble for 
particles is constructed similarly: one imposes the total energy of the system using a J-function on an a priori uniform 
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measure on the canonical variables (position and niorncntuni) . However, in the latter case, this procedure is justified by 
the Liouville theorem which selects the relevant canonical variables. In general, however, there is an ambiguity since the 
assumption of uniformity is not invariant under changes of variables. 

It is instructive to discuss the simplest case whore the Edwards assumption can bo discussed in details, and perhaps tested 
experimentally or numerically. Consider, as proposed by Ertas and Halsey, a single disk in a wedge [64]. In equilibrium, there 
are two contact points and therefore four unknowns: /i,Ar,/i,T and /2,jv,/2,t, where fa,T > means that the force pushes 
upwards. These forces must lead to equilibrium, which gives three equations. There is therefore one degree of freedom which 
is not fixed by the equilibrium requirement, and is dynamically selected. It is easy to see that one must have /i,r = /2,r = /r 
and fi,N = f2,N = Jn- The Edwards measure then reads: 

P{.fN,.fT) = ^s(^fN sint/. + /t costA - ^Mp) e(^/iv - |/t|), (2) 

where ip is (half) the opening angle of the wedge. From this result, one can compute the distribution of the 'mobilization' ratio 
r = /t//jv, which is found to be parabolic: 

P('r) (X (sin-)/) + r cosi/')^ {— H < r < fi) , (3) 

and, of course, zero outside the allowed interval [— /x, fi] . One could test this simple predictions by repeatedly tapping spheres 
made of different materials, and investigate the relevance of the tapping mode and the contact dynamics on the statistical 
ensemble of forces that one generates. This would be quite a valuable starting point, before speculating on more complex 
multi-grain situations. We will discuss below some numerical results that indeed suggest some dependence of the statistics of 
forces on the microscopic dynamics. 

The Edwards prescription is in fact at the heart of the simplest 'scalar' model for the statistics of forces in granular materials, 
which was proposed in [93,46] to account for the empirical exponential tail in the distribution of forces. Although this model 
represents a highly stylized view of granular systems and cannot be expected to be accurate, it provides both an extremely rich 
theoretical benchmark and a pedagogical starting point for more elaborate descriptions. 



II. THE SCALAR MODEL I: DISCRETE VERSION 

A. Definition and motivation 

The drastic simplification of the scalar model is to only retain one component of the stress tensor, namely the 'weight' 
w = azz, and correspondingly, to only consider the force balance equation along the vertical axis. Again for simplicity, one 
can think that the grains reside on the nodes of a two-dimensional lattice, and are labeled by two integers: i in the horizontal 
direction and j in the vertical direction; j increases as one moves downwards. The equilibrium equation can then be written 
as: 

w{i,j) = wo + q+{i- l,j - l)w{i - l,i - 1) + q-{i + 1, j - l)w{i + 1, j - 1) (4) 

where 'wq' is the weight of each grain, and q±{i,j) are 'transmission' coefficients giving the fraction of weight which the 
transmits to its right (resp. left) neighbour immediately below, such that q+{i,j) + q-{i,j) = 1 for all i,j's. The case of an 
ordered pile of identical grains and identical conditions for each contacts corresponds to q± = ^. In this case, the equation 
for the w's become identical to the Master equation describing the population of un-biased random walkers in a one space 
dimension (corresponding to i), evolving in 'time' (corresponding to j): 

w{i,j) = Wo -f ^ [w{i- \,j - 1) +1X1^ + 1,3 - 1)] . (5) 

The term wo is a constant source of particles that arc created uniformly in space and in time. We will explore this analogy 
further below. 

Now, grain assemblies are usually not perfectly ordered: grains have various shapes and sizes; there are packing defects and 
irregularities; even for a perfectly ordered packing of identical grains one can expect that the history has imposed different 
contact loadings. In the above language, it means that the q±{i,j) arc not all identical and reflect the above sources of 
randomness. The idea of Edwards, in this highly simplified framework, can easily be worked out, since each 'blocked' state 
corresponds to a particular choice of - say - q+{i,j) for all Provided q-{i,j) = 1 — q+{i,j), equilibrium is ensured. The 
uniform measure on all blocked states, advocated by Edwards, merely translates, in the present case, as a uniform probability 
distribution for q+ (or q_) between and 1. This defines the q model, which was originally written with an arbitrary number N 
of downward neighbours (A^ = 2 in the example above), and can thus be (in principle) generalized to throe dimensions [93,46]. 
In this case, there are N coefficients qa, a = 1, - ■ ■ ,N per grain, and the Edwards measure corresponds to the choosing all the 
Qa independently on each node i,j such that: 
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(6) 



Therefore, in the present case, Edwards prescription can be exphcitly followed. It may seem a priori that this microcanonical 
assumption is too simple and can only lead to trivial predictions. However, as we discuss below, this is not the case: much 
as for gases where the microcanonical hypothesis can be used to derive, for example, the Maxwell distribution for the particle 
velocities, the g-model predicts a non trivial distribution for the local vertical stress. 



B. Stress distribution and the exponential tail 

The case of a uniform distribution of the g's is interesting because it leads to an exact solution for the local weight distribution 
P{w) for large heights. Let us assume for the moment that the weights on neighbouring sites become asymptotically independent 
( [46,129]). Then P{w) obeys the following mean-field equation: 

Pj+i{w) = / dqidq2p{qi)p{q2) / dwidw2Pj{wi)Pj{w2)6[w - {wiqi + W2q2 + wo)] (7) 
Jo Jo 

where p{q) is the distribution of q, here taken to be p{q) = 1. In the limit j — > oo, the stationary distribution P* of this 
equation can be explicitly constructed and is given by: 

P (w) = ^exp-= (8) 

where 2w = jwo is the average weight. For N ^2, the distribution is instead a Gamma distribution of parameter N; its small 
w behaviour is w'^~^ while the large w tail is exponential. Liu et al. [93,46] have argued that this behaviour is generic and 
survives deviations away from the strict Edwards prescription: for example, the condition for the local weight w to be small 
is that all the A'^ q''s reaching this site are themselves small; the phase space volume for this is proportional to w'^"^ if the 
distribution p{q) is regular around g = 0. However, if instead p{q) oc q''^^ when q is small, one expects P''{w) to behave for 
small w as w^", with a = 1 — Ny < 0. 

Similarly, the exponential tail at large w is sensitive to the behaviour of p(g) around q = 1- In particular, if the max:imum 
value of q is qm < 1, one can study the large w behaviour of P*{w) by taking the Laplace transform of equation (7). One finds 
in that case that P* (w) decays faster that an exponential: 

log P* (w) oc^^oc -w^ P = 7^^^ (9) 

log qKiN 

(Notice that (3=1 whenever qu = 1, and that /J — > (X) when = 1/A'': this last case corresponds to an ordered lattice with 
no fluctuations). In this sense, the exponential tail of P*{w) is not universal: it requires the possibility that one of the q can 
be arbitrarily close to 1. This implies that all other g's originating from that point are close to zero, i.e. that there is a non 
zero probability that one grain is entirely bearing on one of its downward neighbours. 

The success of the g-model with a uniform distribution of the g's is that it provides a simple explanation for the ubiquitous 
exponential tail for the distribution of forces, observed in many experimental and numerical situations. Note in particular that 
this exponential tail was observed in a regular packing of grains, suggesting that very strong heterogeneities in fact exist at 
the contact level. On the other hand, the probability to observe very small w is much underestimated by equation (8): see 
[29,114,13,41]. This might be due to the fact that arching effects are absent in this scalar model. A generalization of the 
g-model allowing for arching was suggested in [40], which dynamically generates some sites where g+ = 1 and g_ = (or vice 
versa). This indeed leads to much higher probability density for small weights. For a more detailed discussion of the relation 
between the g-model the experimental situations, in particular the role of boundaries, see [130]. 



C. The 'critical' case 



There is however one special case of particular interest where the results for P{w) are qualitatively different. Suppose that 
g is a random variable that only takes the values or 1 with probability 1/2. This is called the Takayasu model, which is a 
model for directed river networks for example: at each site of the lattice a river flowing 'south-east' or 'south-west' is randomly 
deflected to the left or to the right. Rivers coalesce upon meeting. The 'source' term wo here describes a constant density of 
'springs' that feed the river network. It can also be seen as a model of diffusing and aggregating clusters in a solution with a 
constant density of 'monomers' (that again play the role of the source term). In this case, it turns out [136] that the stationary 
distribution P*{w) becomes, for large heights, a power law, P*{w) oc w~^~'^, with p, = 1/3 in dimension d = 2 (i.e. one 'spatial' 
and one 'temporal' dimension). 
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The exponent /i was derived analytically but can also understood as follows. Since the direction of the 'rivers' is, at each step, 
a random variable, the typical 'basin of attraction' of a given site is a parabolic object of height t and width y/i. Therefore, on 
the order of t^^^ 'springs' at most can contribute to the river flux on a given site; in other words, one expects the distribution 
P* (w) to take the scaling form: 

where F is a certain function which falls off fast for large arguments. On the other hand, since one must have (w) = wot 
exactly, the exponent fj, is fixed: 

r dwwP*{w)=wot = t^''^'' r dw^M^^^i. (11) 

Jo Jo '^^ 3 

(Note that the integral over u is convergent for this value of yu.) 

Therefore, this model generates a power-law distribution for the local masses, and was proposed early on as a model of 

'self-organized criticality'. In fact, the model is critical in the sense that any deviation of p{q) from a sum of two equal delta 
peaks at and 1 leads to an exponential truncation of P*{w) at large w's. Let us add several remarks: 

• One can compute higher moments of the local weight, to find (w^) ~ ^89/2-1/2 ^ ^ -j^^g particular, (w^) ~ t^^^, a 
result that we will recover below using direct method. 

• One can also generalize the model to higher (spatial) dimensions, where one finds n = 1/2 for all d > 3 (with logarithmic 
corrections in d = 3) . 

• Consider a rectangular sample of width W and height H. What is the order of magnitude of the largest weight encountered 

at the bottom ? For a pure power law distribution such as Eq. (10) with t 00, the maximum value of w is known to be 
of order uimax ^ W^^^ . This estimate can obviously only be valid if Wmax is found to be much smaller than the truncation 
imposed by the function F , which is of order H^^^ . This requires W < ~ H^/'^. However, in this regime where 

W is smaller than the diffusion length, the very argument leading to Eq. (10) breaks down, since the maximum weight 
now scales like WH, and not as H^^^ . Extending the argument, we now find that the distribution of weights reads: 

Therefore, (a) for W > 7?^/^ the maximum value of w is imposed by the cut-off and of order H^^^ , and not A'''^. 
Correspondingly, the participation ratio Y2 = 'Ylfi=\ Wi / {woHW}^ that characterizes the 'localization' of the weight is of 
order H^/'^/W < 1. In the case (b) W <C H^^"^, on the other hand, the weight is localized on a finite number of sites, 
and Y2 ~ 1. 



III. THE SCALAR MODEL II: CONTINUOUS LIMIT AND PERTURBATION THEORY 

A. Continuous limit of the scalar model 

Let us focus on the case N = 2 and define v to be such that q±{i,j) = (1 ± v{i,j))/2. If v is small, the local weight is 
smoothly varying, and the discrete equation (4) can then be written in the following differential form: 

dtw + dx{vw) = p + DodxxW (13) 

where x = ia and t = jr are the horizontal and (downwards) vertical variables corresponding to indices i and j, and a and 

T are of the order of the size of the grains. The vertical coordinate has been called t for its obvious analogy with time in a 
diffusion problem, p is the density of the material (the gravity g is taken to be equal to 1), and Do a 'diffusion' constant, which 
depends on the geometry of the lattice on which the discrete model has been defined. This diffusion constant is of the order of 
magnitude of the size of the grains, o. 

In this model and in the following, we shall assume that the density p is not fiuctuating. Density fiuctuations could be 
easily included; it is however easy to understand that the resulting relative fluctuations of the weight at the bottom of the pile 
decrease with the height of the pile H as H^^^^ , and are thus much smaller than those induced by the randomly fluctuating 
direction of propagation, encoded by q (or v), which remain of order 1 as. H ^ 00. 

Two interesting quantities to compute are the average 'response' G{x,t\xo,to) to a small density change at point {xo,to)i 
measured at point {x,t), and the correlation function of the force field, C{x,t,x' ,t') = {w{x,t)w{x' ,t'))^, where the averaging 
is taken over the realizations of the noise v(x,t). 

Equation (13) shows that the scalar model of stress propagation is identical to that describing tracer diffusion in a (time 
dependent) flow v{x,t). This problem has been the subject of many recent works in the context of turbulence [38]; interesting 
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qualitative analogies with that field can be made. In particular, 'intermittent' bunching of the tracer field correspond in 
the present context to patches of large stresses, which may induce anomalous scaling for higher moments of the stress field 
correlation function. 



• Fourier transforms. 



The limit where a is ill defined and leads to a divergence of the perturbation theory for large wave-vectors k. We thus 
choose to regularize the problem by working within the first Brillouin zone, i.e., we keep all wave vector components within the 
interval I = [—A, +A], where A = ^. Our Fourier conventions for a given quantity / will then be the following: 

f{^,t)=r^e^'yik,t) (14) 

J —A 

+ 00 

f{k,t) = £, e-'*"/(x,i) (15) 

x= — oo 

One has to be particularly careful when computing convolution integrals, such as J ^fi{q)h{k — q) which must be understood 
with limits —A + fe, A (resp. —A, A + fe) if fc > (resp. fc < 0). An important example, which will appear in the response 
function calculations, is: 



Let us then take the Fourier transform of equation (13) along x, to obtain: 



{dt + Dok'')Wk = Pk + ik J ^WqVk-q (17) 



Our aim is to calculate, in the small-fc limit, the average response (or Green) function G{k, t — t') defined as the expectation 
value of the functional derivative {5w{k,t) /5p{k,t')); and the two points correlation function of w, C{k,t) = (w{k,t)w{—k,t)} . 



• The noiseless Green function. 



The noiseless (bare) Green function (or 'propagator') Go is the solution of the equation where the 'velocity' components Vq 
are identically zero: {dt + Dok^)Go{k, t — t')= 5{t — t') which is: 

Go{k, t-t')= e{t - t')e-°o'='(*-*') (18) 
In real space, the propagator is simply the heat kernel, 

n^{^. t-t')= — P~4Do(t-t') (19) 

-^47rDo(t-t') 

This shows that in the non-disordered scalar model, the stress 'difFuses', as already noticed above in the discrete formulation 
of the model, see Eq. (5). 

• Statistics of the noise v{x,t). 

The noise term v represents the effect of local heterogeneities in the granular packing. The mean value of the noise v is taken to 

be zero, and its correlation function is chosen for simplicity to be of the factorable form {v{x,t)v{x' ,t')) — gx{x — x')gt{t — t'), 
where Qx and gt are noise correlation functions along x and t axis. We shall take g^ and gt to be short-ranged (although this 
may not be justified: fiuctuations in the micro-structure of granular media may turn out to be long-ranged due to e.g. the 
presence of long stress paths or arches), with correlation lengths and it. Our aim is to describe the system at a scale L 
much larger than both the lattice and the correlation lengths: a, r, ix, it ^ L. This will allow us to look for solutions in the 
regime k,E ^ 0, where k and E arc the conjugate variables for x and t respectively, in Fourier-Laplace space. However, we 
shall see below that the limit a,T,lx,lt can be tricky, and nmst be treated with care: this is because the noise appears in 
a multiplicative manner in equation (13). For computational purposes, we shall often implicitly assume that the probability 
distribution of v is Gaussian; this might however introduce artifacts which we try to discuss. 

• Ambiguities due to multiplicative noise. Ito vs. Stratonovitch. 



7 



In equation (17), we have omitted to specify the dependence on the variable t. There is actually an ambiguity in the product 
term WqVk-q- In the discrete g-model model [93], the q±'s emitted from a given site are independent of the value of the weight 
on that site. In the continuum limit, this corresponds to choosing Wq{t) to be independent of Vk-q{t), or else that the w's must 
be thought of as slightly posterior to the w's (i.e the product is read as Wq{t — 0)vk-q{t + 0)). In this case, the average of 
equation (17) is trivial and coincides with the noiseless limit; hence G = Gq. This can be understood directly on the discrete 
model by noticing that the Green function G{i,j\0,0) can be expressed as a sum over paths, all starting at site (0,0), and 
ending at site 

G(i,j|0,0)= J2 n ^i^'^'^) (20) 

paths V {k,l)ev 

where the q±{k,l) are either q+{k,l) or q-{k,l), depending on the path. Since each bond q±{k,l) appears only once in the 
product, the averaging over q is trivial and leads to: 

G{i,j\0,0)= 2-^=Go(i,j|0,0) (21) 

paths V 

(Note that this argument fails for the computation of the correlation function C, since paths can 'interfere'. We shall return 
later to this calculation.) 

The above choice corresponds to Ito's prescription in stochastic calculus. Another choice (i.e. Stratonovitch's prescription) 

is however possible, which corresponds to the proper continuum time limit in the case where the correlation length £t is very 
small, but not smaller than a. In this case, the w's and the v's cannot be taken to be independent. This is the choice that we 
shall make in the following. 



B. Calculation of the averaged response and correlation functions. 

Two approaches will be presented. The first one, based on Novikov's theorem, leads to exact differential equations for G 
and C, which can be fully solved. The second one is a mode-coupling approximation (mca), based on a re-summation of 
perturbation theory. It happens that, for this particular model where the noise is Gaussian and short range correlated in time, 
both approaches give the same results, because perturbation theory is trivial. In other cases, though, where exact solutions are 
no longer available, the mca is in general very useful to obtain non perturbative results. 

We shall see that the effect of the noise is to widen the diffusion peak: Do is renormalized by an additional term proportional 
to the variance of the noise v. 



{w{k,t)v{k',t)) = J*dt' J dg^^jg^^ {v{q,t')v{k',t)) 



• Novikov's theorem. Exact equations for G. 

Novikov's theorem provides the following identity, valid if the v arc Gaussian random variables: 

(22) 

Such a term actually appears in equation (17), after transformation into an equation for G: 

(dt + Dok^)G{k, t-t') = 6{t - t') - ^kJ^^ J ^ {v{q, t)w{k - q, t)} (23) 

In the limit where £x = a —^ 0, the noise correlation is of the form: {v{q, t)v{q\t')) = 2iTa^S{q + q')gt(t — t'), with gt peaked in 
t = t' such that f{t')gt{t — t') ~ f{t)gt{t — t') for any function /. Prom formally integrating equation (17) between t' and t, 

one can express the equal-time derivative Sw/Sv as: 



Sw{k, t) 



= -ikw{k - k' ,t) (24) 

t'=t-o 



Sv{k',t') 

and thus obtain: 

dq 



{dt + Dok^)G{k, t-t') = 5(t - t') - 2wa^kG{k, ^ " ^ d'^'dti^' " J 



2^(fc-9) (25) 



Using the shape of the function gt, the first integral is 1/2. The second one is a convolution integral, and its value is 
Ak/2n + 0{k^) (see equation (16)). The final differential equation for G is then, in the small-fe limit, a diffusion equation with 
a renormalized diffusion constant: 
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Dr = Do + — (26) 

It is interesting to note that the model remains well defined in the limit where the 'bare' diffusion constant is zero, since a non 
zero diffusion constant is induced by the fluctuating velocity v. This would not be true if Eq. (13) was interpreted with the 
Ito convention, where the fluctuating velocity would not lead to any spreading of the average density. 

The most important conclusion is thus that, in the present scalar model, stresses propagates essentially vertically, even in 
the presence of disorder: taking I ^ a (where tx ^ (-t ^ (), the response at depth H io a. small perturbation is confined within 
a distance oc ^JDrH from the vertical. Since Dr ~ £^/a, ^JDrH is much less than H in the limit where H ^ £^/a, i.e. when 
the height of the assembly of grains is much larger than the grain size. 

• Exact equations for C. 

Exact equations can also be derived for C, following very similar calculations. From equation (17), one can deduce 
the corresponding one for w{k,t)w{—k,t). Upon averaging, Novikov's theorem has to be used on quantities such as 
{w{k,t)v{q,t)w{—k — q,t)), finally leading to: 



{dt + 2Diik^)Cik,t)=a^k^ 

where = (27r)cr^. Going to Laplace transforms in time leads to: 

(E + 2DRk'^)C(k, E) - C{k, t = 0) = B^k" [C{E) + 2p^E-^] (28) 

where C{E) = J ^C{k,E) = J dte~'^''C{x = 0,t). This allows one to get a closed equation on C{E): 

CiE) = C{k,t = 0) fdk a'k' ,g , 3^ _ (29^ 

^ ' J 2'kE + 2DRk^ J 2-n:E + 2DRk^ \ ^ ' ^ ) ^ ' 

In the limit where — > 0, the first term on the right hand side is of order E~^^'^, whereas the coefficient of C{E) is equal to 
/2Dr. This shows that one has to distinguish two cases: 

• < 2Dr. This corresponds, for the discrete model, to the generic case where p{q) is not the sum of two delta peaks at 
9 = and 9 = 1. The equation for C{E) becomes, for _E — > 0: 

- -TT- ) C(E) « -^2p2£;-^ (30) 
2Dr J ^ ' 2Dr ^ ' ^ > 

or C{E) ~ E~^ . Transforming back to real times leads to C(a; = Q,t) ~ , which is consistent with the results obtained 
within the discrete scalar model, where the local weight is of order i, with relative fluctuations of order one. Once one 
knows C{x = Q,t), one can obtain, from Eq. (27), the full function C{x,t). When a; S> a, one finds: 

C{x,t)=t'''c(^^, (31) 

where C(.) is a scaling function computed explicitly in [115]. This results shows (i) that the correlations extend over 
distances \/t, as expected from the diffusive nature of the model, and (ii) that C{x ^ 0,t) <^ C(x = 0), meaning 
that asymptotically correlations between different sites vanish. The assumption that different sites become independent 
(which is a stronger statement) was already used above to obtain the Gamma distribution of the local weights in the 
scalar model. We see that this assumption is at least consistent. 

cP' = 2Dr. This corresponds exactly, as shown in [115] to the critical Takayasu model where q can only take the values 
and 1 with equal probability. In this case, a more careful analysis of the limit £ — > must be performed. Using the 
fact that: 



gc(g,t) + pV 



(27) 



f dk 2DRe _ fdk 1 r-j— 

J 2^E + 2DRk^-^-^J 2;^i; + 2D«fc2 -l-V-B/81)fl, 



(32) 



we now find that C{E) ~ £ "^^^j ot C{x = 0,t) ~ t^^^. Again, this is consistent with the direct scaling analysis presented 
above. Extending the analysis to x ^ now leads to: 

C{xT^O,t)=fcJ^], (33) 



where the critical scaling function Cc(.) can also be computed explicitly [115], and is different from C(.). Note that the 
asymptotic de-correlation of different sites still takes place, since C{x ^ 0,t) <^ C{x = 0,t). 
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• Perturbation theory. 

The above method gives exact results, essentially because v{x,t) is short range correlated in time: 6w/Sv is then only needed 
at coinciding times, where it is known, and equal to 1. This would not be true in general; furthermore Novikov's theorem 
requires v to be Gaussian. It is thus interesting to show how a systematic perturbation scheme can be made to work by the use 
of diagrams to represent equation (17). The MCA (Mode Coupling Approximation) is then a particular re-summation scheme 
of this set of diagrams, which was discussed in detail in [19], which sometimes provide interesting non perturbative results. 

Equation (17) is multiplied on the left by the operator Go (see equation (18)), and then re-expressed as follows: 

w{k,t) = Go{k,t) ® p{k,t) - ikGo{k,t) ® j ^w{q,t)v{k - q,t) (34) 

(gi meaning a f-convolution product. This equation can be represented with diagrams as follows: as shown in figure 1, we 
represent the source p by a cross, the 'bare' propagator Go by a plain line and the noise v by a dashed line. The first term of 



p{k, t) = x Go{k, t) = — L— wo{k, t) = L-x 

v(k^t) = ---i-- G(k,t) = * w(k,t) = 

k — q ^''^ q — k 



t' 



FIG. 1. definition of various diagrams. 



equation (34), which is the noiseless solution wo, is then obtained as the juxtaposition of a plain line and a cross. The arrow 
flows against time (i.c it is directed from t to t' < t). The juxtaposition of two objects means a f-convolution product. By 
definition w is represented by the juxtaposition of a bold line and a cross (this is consistent with the identification of a bold 
line with the full propagator G). The diagrammatic version of equation (34) is then: 

A k-q 

I 

k , , k ^ , k ^ q ^ , 
* — K = > — X + > — ' >*-^ 



(35) 

The 'vertex' stands for —ik J the two emerging wave vectors being q and k — q (node law). One can now iterate this 
equation. To second order, one obtains: 



* — K = + * ' > X + ^ ' — ^X (36) 

The corresponding equation for G is obtained by taking the derivative S/5p, and averaging over the noise v. Since {v} = 0, the 
second diagram vanishes. We represent the noise correlator by a dashed line with a centered circle (see figure 1), and obtain: 



(37) 



or G = Go + GoSGo, where E is called the self-energy (see figure 1). Actually, one can re-sum exactly all the diagrams 
corresponding to GoSGo, GqSGoSGo to obtain the Dyson equation G = Go -|- GoSG. 

The MCA amounts to replacing the 'bare' propagator in the diagram for E by the full propagator G. (Note that the MCA is 
of course exact to second order in perturbation theory). We then obtain a self-consistent equation for G: 

EmCA = Gq — GmcA = ' * (38) 
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FIG. 2. Example of a diagram included in the mca. Note that this diagram is zero if the noise is (5-correlated in time. 

Diagrams like the one drawn in figure 2 are now also included. The self-energy Emca can be easily computed, we get 

E„cA(A;,t - t') = -2na^k J ^qGucK{q,t - t')gt{t - t') (39) 

In the special case where gt is peaked around t = t' , we can make the approximation G{q, t — t')gt{t — t') ~ G{q, Q)gt{t — t') = 
gt{t — t') (since by definition G{q,Q) = 1). We thus get, using equation (16)), Emca(/c, t — t') = —a^Ak^gt{t — t'). The 
expression for Gmca is thus identical to the one obtained with the exact approach, as can be seen by comparing equation (25) 

and Go = 1 + EG. 

Note that one can also calculate the influence of a non zero skewness ?, or kurtosis n, which are the normalized third and 
fourth cumulant of the noise v. In this case, three and four dashed lines (corresponding to v) can be merged, leading to a 
contribution to D, of the order of ?cr^ or /tcr*. 

Let us turn now to the calculation of the correlation function C{k, t) = {w{k, t)w{—k, t)}. The basic object which corresponds 
to the self-energy is now the 'renormalized source' spectrum S{k, t, t') defined as: C = G® S ®G. S is drawn as a filled square. 
Sq (empty square) is the correlation function source term which encodes the initial conditions (see below). The two first terms 
of the expansion are 

/ \ 

/ \ 

-m- = ^ + A h +... 

^-or (40) 

Here again, we transform the perturbative expansion into a closed self-consistent equation for S by replacing Go and So in (40) 
by G and S respectively. The final equation for G reads: 



/ 

/ 

(41) 



or, written explicitly. 



C{k,t)= I dt' I dt"G{k,t-t')So{k,t',t")G{-k,t-t")+ 
Jo Jo 



a^fc^^ dt' dt"G{k,t-t') J ^C{q,t'X)at{t' -t")G{-k,t-t") (42) 

If we choose the source term to be an overload localized at t = 0, we get: So = {p{k,t')p{—k,t")) = C{k,0)5{t')5{t"). 

Using the fact that gt is peaked around t' = t" , we again recover exactly the equation (27) above, showing again that mca 
is exact in this special case. 

C. Further results: the un-averaged response function 

The average Green function described above is thus a Gaussian of zero mean, and of width growing as s/Dnt. However, for 
a given environment, the Green function is not Gaussian, presenting sample dependent peaks (see Figure 3). Note however 
that, contrarily to what we shall find below for the tensorial case, the un-averaged Green function remains everywhere positive. 
Furthermore, the quantity [a;](t), defined as the displacement of the centroid of the weight distribution beneath a point source 
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FIG. 3. Averaged (bold lino) and un-avoraged (thin lino) rosponso function of the scalar model, obtained numerically by 
simulating the q-model. One can notice how 'non self averaging' is the response function, i.e. how different it is for a given 
environment as compared to the average. Note also that the un-averaged Green function is everywhere positive. 



in a given realization: 

l*) = />'«'^ («) 

typically grows with t. More precisely, one can show that: 

{[x]{t)) = but {[xf{t))o:t^^^ (44) 

meaning that the 'center' of the Green function wanders away from the origin in a sub-diffusive fashion, as t^^*. This behaviour 
has actually been obtained in another context, that of a quantum particle interacting with a time dependent random envi- 
ronment. Physically, the g-model can indeed be seen as a collection of time dependent scatterers, converting in-going waves 
into outgoing waves with certain partition factors g+ , q_ = 1 — g+ , a problem equivalent to the one dimensional Schroedinger 
equation with a time dependent random potential (sec the discussion in [122,47]). In two dimensions (plus time), the wandering 
of the packet center [x]{t) is only logarithmic (and disappears in higher dimensions [47]). 

Similarly, the participation ratio Y2 = J dx (w'^{x,t)'^ can be computed, and is found to be ~ t~^^'^, which means that the 
weight is not localized on a finite number of sites in this model when t —>■ 00. 



D. The scalar model with bias: Edwards' picture of arches 

Up to now, wc have considered the mean value of v to be zero, which reflects the fact that there is no preferred direction for 
stress propagation. In some cases however, this may not be true. Consider for example a sandpile built from a point source: the 
history of the grains will certainly in-print a certain oriented 'texture' to the contact network, which can be modeled, within 
the present scalar model, as a non zero value of (v), the sign of which depends on which side of the pile is chosen. In other 
words, the isotropic Edwards assumption for the local stress transmission is expected to break down when the history of the 
packing explicitly breaks a symmetry. 

Let us call Vb the average value of v on the a; > side of the pile, and — Vb on the other side. The differential equation 
describing propagation now reads, in the absence of disorder: 

dtw + dx [Va sign(x)w] = p + Dod^xW (45) 

For a constant density p = po, and for Do = 0, the weight distribution is then the following: 
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w{x,t)= ^ forO<a;<Vbi 

w{x,t) = iorVot<x<ct (46) 

where c — 1/ tan</i {(f) is the angle made by the slope of the pile with the horizontal x axis). For Do 7^ 0, the above solution is 
smoothed. In any case, the local weight reaches a minimum around a; = 0. Equation (45) gives a precise mathematical content 
to Edwards' model of arching in sand-piles [59], as the physical mechanism leading to a 'dip' in the pressure distribution [128]. 
As discussed elsewhere [145,146], this can be taJsen much further within a tensorial framework (see Section V A). 

Equation (45) with noise can in fact be obtained naturally within an extended q'-model, with an extra rule accounting for 
the fact that a grain can slide and lose contact with one of its two downward neighbour when the shear stress is too large [40] . 
This generically leads to arching; in the sandpile geometry and for above a certain probability of (local) sliding, the effective 
'velocity' Vb becomes non zero and the weight profile (46) is recovered [40]. However, this extra sliding rule implicitly refers to 
the existence of shear stresses, which are absent in the scalar model, but which are crucial to obtain symmetry breaking effects 
modeled by a non zero Vb (see also the discussion in section VA. It is thus important to consider from the start the fact that 
stress has a tensorial, rather than scalar, nature. This is what we investigate in the following sections. 



IV. STATIC INDETERMINACY; ELASTICITY AND ISOSTATICITY 

A. Elasticity and response functions 

As mentioned in the introduction, the sole equilibrium equations are not sufficient to determine the stress tensor of an 
arbitrary material. In d = 2, one has two equations and three independent components of the stress tensor. In d = 3, there 

are three equations for six independent components of aij. In clastic materials, this indeterminacy is lifted when one adds the 
constraint that the stress tensor is linearly related to the strain. The most general linear relation between stresses and strains 
is given by: 

Uij = ^ijklUkl (47) 

where aij denotes the components of the stress tensor, m is the displacement field, Uij = ^{djUi + diUj) those of the strain 
tensor, and summation over repeated indices is implied. The four index tensor \ijki satisfies certain symmetry conditions [89]. 
In order to close the problem for the stress tensor, one imposes a condition of 'compatibility', which in d = 2 reads: 

dluxx + dlu^z - 2dxdzU:rz = 0, (48) 

resulting simply from the fact that the tensor Uij is built with the derivatives of a vector m. This is enough to find a closed 
equation for the stresses (in d = 2) [110]: 

{dt + tdt + 2rdldl) aij = (49) 

where the two independent coefficients t and r can be expressed in terms of the components of Xijki- Isotropic elasticity 

corresponds to r — t = 1. 

Expanding the stresses in Fourier modes, it is easy to see that the solutions of the equations (49) are of the form 



crtt 



dQ^a,(g)e'«^+^^'=«^ (50) 

k 



= C,z- dg^a;,(g)X,e'«^+'^''''^ (51) 

/• + OC 

= a.+ / d(z^afe(g)Xfc'e'^"+'^'=«^ (52) 



where Cxx and Cxz are constants. Prom equation (49) we see that the Xk are the roots of the following quartic equation 

+ 2rX^ + t = 0, (53) 

which has four solutions: 



X = ±\l ~r±y/r^ -t. (54) 

Hence the index k runs from 1 to 4. The four functions ak{q) and the constants Cxx and Cxz must be determined by the 
boundary conditions. A particularly interesting boundary condition is when one imposes a localized force at the top surface of 
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the material. The shape of the stress response function to such a locahzcd force will be of central importance in the following 
discussion. One can establish the existence of various 'phases' in the r, t plane in terms of the shape of the response function, as 
obtained from the calculations presented in [110]. In that plane, the line t = r^, for r < 0, separates the so called hyperbolic and 
the elliptic regions. For t > r^, the above roots are complex whereas for t < r^ and r > the roots Xk are purely imaginary. 
These two regions correspond to the elliptic regime, which is in fact the only accessible one in the context of classical elasticity 
where the coefficients r and t are constrained by the fact that the undeformcd state is a minimum of the clastic energy. As 
shown in details in [110], the response function has a unique, broad peak of width growing linearly with depth in the region 
t < r^ and r > 0, whereas the response function becomes double peaked in the region r < 0, t > r^ (with a width again scaling 
linearly with depth). As one approaches the line t = r^, the two peaks become narrower and narrower before becoming two 
deltarfunction peaks exactly on the transition line. At this point and below the transition, the system is hyperbolic; this limit 
behaviour will actually emerge naturally below in the context of granular materials. 

B. Indeterminacy at the grain level and isostaticity 

Elasticity theory can also be seen as the long-wavelength description of a network of beads and springs, for which the local 
equilibrium equations are fully determined at each node. When the network is disordered, the theoretical difficulty is to compute 
the effective elastic constants in terms of the probability distribution and correlations of the microscopic springs. The same 
problem arises when one wants to compute the effective conductivity, or the effective permittivity, etc. of a composite material. 
But in all these problems, the microscopic equations are sufficient to solve the problem in principle. 

For an assembly of grains, this is not the case. The indeterminacy of the static equilibrium exists already at the grain level 
(see the simple case discussed in section IB 2). In principle, one should describe in details the microscopic history of each 
contact in order to determine the precise configuration of forces within a given packing. There arc however special cases where 
this is not the case, and where all contact forces axe fully determined by the packing geometry. These situations are called 
isostatic, and play a special role. These equilibria are in some sense critical since the opening of one contact necessarily leads 
to some rearrangements. Some arguments have been put forward to suggest that an assembly of grains relaxing towards static 
equilibrium will most probably stop as soon as they are stable, i.e., in one of these isostatic states [14]. A similar statement 
is actually exact in the context of mean-field spin-glasses, where the equilibrium states reached dynamically are marginally 
stable, in the sense that the spectrum of the eigenvalues A of the Hessian (matrix of second derivatives of the energy) vanishes 
precisely at A = 0. Before discussing the validity of the idea that these marginal states play a special role in granular media, 
let us first discuss the geometrical conditions necessary for isostaticity [125]. 

Consider frictionless grains in two dimensions. There are, per grain, two equilibrium equations since the torque is auto- 
matically zero, and one (normal) force per contact that must be determined. If is the total number of grains and Nc the 
total rmmbcr of contacts, the number of unknowns is Nc and the number of equations in 2A'^. Therefore, one can (generically) 
find static solutions only if Nc > 2N. Since each contact concerns two grains, the average number of contact per grain is 
Tic = 2Nc/N and the condition for the existence of solutions is ric > 4. The marginal case is when ric = 4, where the number 
of unknowns is equal to the number of equations, and corresponds to the isostatic case. The hyper-static case corresponds to a 
strict inequality. If the friction coefficient is non zero, then the zero-torque condition provides a third equation for each grain. 
If we call ip the fraction of contacts where friction is fully mobilized (i.e. such that \Ft\ = (j-Fn), one has (p ■ Nc + (1 — ip) ■ 2Nc 
unknowns (one per mobilized contact, two for un-mobilized contacts). The stability condition now reads nc(l — (p/2) > 3. 
(Note that the frictionless case corresponds to tp = 1). In three dimensions, the same argument leads to (3 — ip)nc > 12, the 
isostatic case corresponding to an equality. The corresponding 'stability' diagram is shown in Fig. 4. 

C. Numerical simulations and Edwards' assumption 

So, are three dimensional packings of grains, obtained by letting the grains lose their kinetic energy and come to rest, 

generically isostatic? The only quantitative study we are aware of is that of Silbert et al. [125], where these authors perform 
Molecular Dynamics simulation of monodisperse grains with a specific form of contact dynamics and a certain energy dissipation 
coefficient at each collision. The results reported in [125] are compatible with isostaticity for frictionless spheres, /i = 0. In this 
case, the equilibrium packings are indeed found to obey the condition ric = 6, as expected from the general results of [118,102]. 
However, when the friction coefficient is non zero, these authors find that the static configurations are such that (a) the fraction 
of fully mobilized contacts is = and (b) the rmmbcr of contacts per grain seems to saturate, in the limit of hard grains, to 
a value of tIc > 4, suggesting that the packing is not isostatic. The value of Uc appears to depend significantly on the value of 
the friction coefficient and the restitution coefficient; it appears from their data that smaller restitution coefficients (i.e. more 
damping) lowers the value of ric, perhaps down to the isostatic limit for large damping. This dependence on the details of the 
dynamics indicates that no universal statement about the statistics of 'native' packings (i.e. obtained without further tapping) 
can be made. Again, the simplistic situation discussed in section IB 2 would provide a useful benchmark. 

What happens when one of these native states (possibly hyperstatic) is vibrated? Does the packing wander inside the stable 
regime (see Fig. 4) or remains near the isostatic boundary ? Here again, it is useful to recall the results of mean field p-spin 
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FIG. 4. Stability diagram in the plane 99, ric for three dimensional assemblies. The plain line is the isostatic line that separates 
stable packings from unstable packings. The point <p = 1, ric = 6 corresponds to frictionless particles. The cross for = 
represents the numerical result of [125]; the dotted arrow is a possible path of the packing when vibrated. 



glasses, where the 'vibrations' (temperature) keep the system along the ridge of marginally stable states. In this case, the 
reason is the exponential dominance of the number of these states over the 'deeper' (more stable) ones. Therefore, even if 
blocked states are a priori equiprobable (as postulated by Edwards), the most probable situation is to observe the system in a 
marginal state. If this argument can be transposed to granular packings, then ideas that Edwards expressed in different contexts 
(i.e; that blocked states are equiprobable and that only marginal (isostatic) states are important [14]) would be reconciled. It 
would be very interesting to compute the number of metastable states as a function of the isostatic index ric in some (possibly 
artificial) model (on this point, see the attempt in [99]). 



V. A STRESS-ONLY APPROACH TO GRANULAR MEDIA 



We now turn to the discussion of some plausible 'stress-only' closure schemes for the static equilibrium of granular materials. 
We first start by a natural generalization of the scalar g-model to account for the vectorial nature of the forces. Then we show 
how the results of this 'vectorial' g-model can be interpreted more generally in terms of symmetry arguments. 



A. A vectorial g-model 

It is useful to start with a simple toy model for stress propagation, which is the analogue of the scalar model presented in 
section II. We now consider the case of three downward neighbours (see figure 5), for a reason which will become clear below. 
Each grain transmits to its downward neighbours not one, but two force components: one along the vertical axis t and one 
along X, which wo call respectively Ft{i,j) and Fx{i,j)- For simplicity, we assume that each 'leg' emerging from a given grain 
can only transport the force parallel to itself. This assumes that each contact is frictionless. More general transfer rules can be 
considered, where the forces axe chosen, d la Edwards, randomly within the space of solutions - see e.g. [63,131,105,28]. 

We assume that a fraction p of the vertical force travels through the middle leg. Correspondingly, a fraction g+ = (1 — p)(l-|- 
e)/2 travels through the right leg, and a fraction g_ = (1 — p)(l — e)/2 through the left log, so that p + g+ + = 1. (Some 
anisotropic generalizations will be discussed below). Since the total force in these legs is oriented in the direction ip, where is 
the angle between grains, defined in figure 5, the balance of the horizontal force on the grain i, j imposes a specific value for e: 

^- ' ^^(■''^■) (55) 



(1 -p)tanV F,{i,j) 



We therefore discover a very important consequence of the existence of a shear stress F^, which was totally absent from the 
scalar g-model or had to be put by hand: any non zero value of necessarily leads to q+ ^ q- and biases the propagation of 
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FIG. 5. Three neighbour configuration. Each grain transmits two force components to its downward neighbours. A fraction 
p of the vertical component is transmitted through the middle leg. 



Fz in the direction of Fx- This coupling between the two components of the force is at the origin of all the interesting physics 
discussed below. 

Using the value of e to compute the propagation of the forces from one layer to the next, we find: 
F4i,j) = \ - 1, j - 1) + Fx{i + 1, j - 1)] 

+^(1 -rttan^A - 1, j - 1) - Ft(i + 1, j - 1)] (56) 

Ft{i,j) = wo +pFt(2,j - 1) + i(l - p) [Ft{i - l,j - 1) + Ft{i + 1, j - 1)] 

Taking the continuum limit of the above equations leads to: 

dtFt + d^Fx =p + dlFt (58) 

dtF.+d.[clFt] =^dLF. (59) 

where Cq = (1 — p) tan^ Tp, a is the size of the grains, and we have kept the second order diffusion terms, which would be the only 
remaining term in the isotropic scalar description, and plays the role of a smoothing term for the (singular) solutions found 
below. Comparing Eq.(58) with Eq. (45) (with w = Ft), we see that the 'velocity' term introduced by hand in the Edwards 
model is indeed a consequence of the local shear. 

Eliminating (say) Fx between the above two equations loads to a wave equation for Ft (up to a diffusion term which becomes 
small in the large scale limit), where the vertical coordinate t plays the role of time and co is the equivalent of the 'speed of 
light'. In particular, the stress does not propagate vertically, as it does in the scalar model, but rather along two rays, each at 
a non zero angle ±ip such that co = tani^. Note that ip ^ tp m general (unless p = 0); the angle at which stress propagates has 
nothing to do with the underlying lattice structure and can take any value depending on the local rules for force transmission. 
The three-log model was chosen to illustrate this particular point: the number of legs is irrelevant in the large scale limit, and 
the wave structure of the resulting equation is universal. 



B. A constitutive relation between stress components 

The above derivation can be reformulated in terms of classical continuum mechanics as follows. Considering all stress tensor 
components aij, the equilibrium equation reads, ^ 

dtcrtt + dxCTxt = p (60) 
dtCTtx + dxCTxx = (61) 

Identifying the local average of Ft with att and that of Fx with atx, we see that the above equations (58, 59) are actually 
identical to (60, 61) provided atx = <Jxt (which corresponds to the absence of local torque) and, more importantly, that there 
exists a relation between the vertical and horizontal components of the stress tensor: 

(Jxx = Co ou (62) 
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This relation between normal stresses was postulated in [18] as the simplest "constitutive relation" among stress components, 
obeying the correct symmetries, that one can possibly assume. The term "constitutive relation" normally refers to a relation 
between stress and strain, but the model under discussion has no strain variables defined. This particular choice can be 
interpreted as a local Janssen approximation [82], in analogy with the assumption made by Janssen in 1895 in order to describe 
stresses in silos, where the average vertical stress at a given altitude is postulated to be proportional to the average horizontal 
stress. We return later to a more detailed discussion of this type of closure equations. In the present case, the parameter Cq 
encodes relevant details of the local geometry of the packing (friction, shape of grains, etc.) and may thereby depend on the 
construction history of the grain assembly. Only for simple, 'homogeneous' histories (such as constructing a uniform sand-bed 
using a sieve) will c§ be everywhere constant on the mesoscopic scale. Even then, unless an ordered packing is somehow created, 
local fluctuations of Cq will always be present. The influence of these fluctuations will be analyzed below. 

C. Some simple situations 

Using Eq.(62) in Eqs. (60, 61) we find that these can be rewritten as: 

duicTxt- coatt) = (63) 

dviffxt + cofxt) = |, (64) 

where u = x — Cot and v = x + Cot. This shows that in this framework, the forces are transported along the characteristics, which 
is what the word 'force chain' usually implies [61,35]. More precisely, the solution of Eqs. (60, 61) is obtained by projecting 
any 'initial' force (i.e. present at altitude t = 0) onto the two characteristics, and propagating and augmenting each component 
by an amount {p/2)dt, independently along these characteristics. 

The simplest situation is that of an infinitely wide layer of sand, of depth H , with a localized (J-function) overload at the top. 
The additional weight at the bottom then defines the response function of the wave equation, which, in two dimensions, is the 
sum of two 5 peaks localized at a; = ±coH. These 6 peaks are actually diffusively broadened by the second order terms present 
in Eqs. (58, 59). This two-peak response function is notably different from the response function of an isotropic elastic body, 
for which the response function is a single hump of width proportional to the height H. However, for anisotropic elasticity, a 
two peak response function can be observed [73,68,110], but both peaks have a width scaling linearly with H, and not as H^^'^ 
for an hyperbolic equation. The question of the response function is therefore of crucial importance, and will be discussed in 
details below. Recent attempts to determine the response function experimentally seem to favor an elastic like response, at 
least for the strongly disordered systems (see section VI). 

Next, one can consider the sandpilc geometry. For a pile at repose, the position of the free surfaces are x = icz, where 
c = cot (j) with the repose angle. On these surfaces, all the stresses vanish. This boundary condition is then (for given co and 
c) suflicient to solve for the stress field everywhere in the pile. One then finds that the vertical normal component of the stress 
is piecewise linear as a function of x. In particular, for ~coH < x < cqH, att is constant. Therefore, in two dimensions, this 
model [18] predicts a flat-topped stress profile rather than a hydrodynamic pressure hump or a dip. Such a flat top is in fact 
observed when building a pile from a uniform rain of grains. 

For a pile created by depositing grains from above (for example by sieving sand onto a disc) it is natural to expect the free 
surface to be a slip plane. (This is a plane across which the stress components saturate the Mohr-Coulomb condition - see 
below.) Interestingly, this provides a relation between co and the friction angle (f), which reads: Cq — 1/(1 + 2tan^(/)) (note 
that since c = l/tan<i), one has automatically c > co). Under these conditions one finds that the 'plastic' region (where the 
Mohr- Coulomb condition is saturated) extends (in two dimensions) inward from the surface to encompass the outer 'wings' 
of the pile (i.e. cqz < \x\ < cz). This follows from the solution of the model and is not an a priori assumption, of the kind 
commonly made in elastoplastic modeling (e.g., [33,123]; see also the discussion in [104]). In three dimensions, a second closure 
relation is required [18] , but in all cases the stress profile has a broad maximum at the center of the pile. Now, however, the 
Mohr-Coulomb condition is only saturated in the immediate vicinity of the free surface - the 'plastic' region has zero volume 
in three dimensions [18,146]. 

D. Symmetries and Constitutive Relations 

Although above it was motivated in the context of a specific microscopic model of force transfer, the linear constitutive 
relation (62) can be viewed, independently of any microscopic model, as the simplest closure equation compatible with the 
symmetries of the problem. The latter include a local reflection symmetry in which a; — a;o is changed to xo — a: (with xo an 
arbitrary reflection plane) and also a form of "dilational" symmetry on the stress known as RSF ("radial stress field") scaling. 
RSF scaling depends on the absence of any characteristic stress scale, which follows if the Young's modulus of the grains is 
sufficiently much larger than any stresses arising in the granular assembly being studied. Such scaling, which requires the stress 
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distributions beneath piles of different heights to have the same shape, is quite well confirmed in some (but not all) experiments 
on conical sand-piles [f 28,29,123]. 

Even with these two symmetries, one can consider more complicated (nonlinear) constitutive relations among stresses, which 
must be of the form [18]: 



cTxx = clatt T ( ^ ) . (65) 



Note that in general, such a constitutive relation violates rotational symmetry, and can only describe an anisotropic pile (for 
example, built from a rain of grains under gravity) . The only rotationally invariant stress only closure scheme can only involve 
the two invariants of the stress tensor, namely Tr cr and Det cr. Prom dimensional analysis, this relation can only be of the form: 

cos^ 0(Tr af = 4Det cr, (66) 

where the specific choice of the coefficient is such that we recognize the usual Mohr-Coulomb relation. Eq. (66) indeed means 
that there exists a choice of orthogonal axis m, n such that: 



= tan 4> = fi, (67) 

where is the internal friction coefficient. One can easily show that this relation is of the general form Eq. (65), for a particular 
form of J^, since Eq. (66) can also be rewritten such as: 



Co J^{u) 



sin^ + 1 ± 2 sin — cot^ 



(68) 



Viewed as a constitutive equation, the Mohr-Coulomb relation defines a rigid-plastic model whose physical content is to 
assume that, everywhere in the material, a plane can be found across which slip failure is about to occur, hence the name 
"incipient failure everywhere", (ife), given to this model [18,145,146]. 

All closures of the form (65) lead to hyperbolic equations for stresses, although in the general case the characteristic directions 
of propagation (the 'light rays' of the corresponding wave equation) depend on the loading and therefore vary with position. 

An interesting situation arises when local reflection symmetry is broken. This is the case, for example, in sand-piles created 
by pouring from a point source onto a rough surface - which is the usual mode of construction. In such a pile, all grains arriving 
at the apex of the pile roll (in two dimensions) either to the right or to the left. The two halves of the pile therefore have 
different construction histories that are mirror images of each other. This violates local reflection symmetry, and in general 
permits constitutive equations such as: 

2 ^ f sign{x)a^t\ 
(Txx = Co(Tu Q — (69) 

V / 

The simplest case (found e.g. by expanding Q to first order in the shear to normal stress ratio) corresponds to a family of 
(quasi-) linear constitutive relations [146]: 

cTxx = clcTtt + V sign(a;)(Txt (70) 

The previous, symmetrical, case has v = Q. For nonzero v, (70) again leads to a wave equation, although this time anisotropic, 
in the sense that the two characteristic rays make asymmetric angles to the vertical axis. Note also that a; = is a singular 
line across which the directions of propagation change discontinuously. 

Microscopically, v leads to an unequal sharing of the weight of a grain between the two characteristic rays propagating 
downward from it. Such a model can indeed be obtained from rules such as those in figure 5 simply by having an asymmetric 
rules for partitioning the forces between the supporting grains, for example choosing two different angles V+ saiA The 
symmetric case corresponds to = and v oc tani/;+ — tani/)_. For < 0, most of the weight travel outwards on the leg with 
the smallest opening angle; this provides, within a fully tensorial model, a mathematical description of the tendency to form 
arches, as developed by Edwards for the scalar case. 

Solving these anisotropic wave equations for sand-piles in two dimensions one again finds for att a piecewise linear function, 
which now has a maximum at x = when u > Q, but a minimum for < 0, in accord with the arching picture mentioned 
above. If one furthermore imposes, as above, that the free surfaces are slip planes, one finds a relation between Cq, y and (j). 
One again finds the result that the material throughout the outer wings of the pile (exterior to the triangle formed by the 
characteristics passing through the apex) are at incipient (Mohr-Coulomb) failure. 
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E. Boundary conditions and 'fragility' 



All the above closure schemes lead to hyperbolic equations, which crucially differ from the elliptic equations encountered in 

elasticity theory when boundary conditions are considered. Hyperbolic equations indeed require only 'half of the boundary 
conditions to bo specified, the other half being determined by 'propagating' these known boundary conditions along the char- 
acterisitcs through the system. On the contrary, an elliptic equation (such as Laplace's equation) requires the stress (or the 
displacement) on all the surfaces of the body to be specified. If one insists on applying all the boundary conditions appropriate 
to an elastic body, then in general no solution will exist for the hyperbolic equations that correspond to a particular choice 
of constitutive relation. If those boundary conditions arc 'incompatible' in this sense, then within an hyperbolic model, the 
material ceases to be in static equilibrium. This is not different from the corresponding statement for a fluid; if boundary 
conditions are applied that violate the conditions for static equilibrium, some motion will result. Unlike a fluid, however, for 
a granular medium we expect such motion to be in the form of a finite rearrangement rather than a steady flow. Such a 
rearrangement will change the micro-texture of the material, and thereby alter the constitutive relation among stresses. One 
expects it to do so in such a way that the new (spatially inhomogeneous) constitutive relation becomes compatible with the 
imposed forces. 

Within this modeling approach, a granular assembly is therefore able to support some, but not all, of the surface loads that 
would be supportable by an elastic medium. Such models may therefore provide an interesting paradigm for the behaviour of 
"fragile matter" [35]. 

VI. EXPERIMENTAL AND NUMERICAL DETERMINATION OF THE STRESS RESPONSE 

FUNCTION 

The stress 'response function' to a localized overload is of prime interest both from a fundamental point of view but also 
for many engineering applications. On large scales, the extra stresses created by a house within the soil beneath it are indeed 
related to this response function. Therefore, this problem has received considerable theoretical attention in the engineering 
community, where, as mentioned in the introduction, the granular material is often assumed to be a (possibly anisotropic) 
elastic material. Note that an elliptic response function corresponds to a favorable case for stability since the stresses are 
efficiently dispersed in space, whereas a hyperbolic response function would lead to a rather localized stress field. 

In spite of its importance, the response function of granular assemblies has only very recently been measured experimentally 
[69,116,103]. Various methods were used: one is a direct quantitative measurement of the response at the bottom of a 3D 
packing using a local stress probe based on the deformation of a hard membrane [116], another method is based on carbon 
paper imprints created by a monodisperse 3D packing [103]. For 2D packing the photo-elastic response of polymeric grains 
was used [69] in order to evidence the inter-particle force path. This is a semi-quantitative method but it allows to visualize 
directly the response in the bulk as well as the topology of the path followed by the stress chains. Such an observation will 
be used to built the theoretical proposition exposed in the following section. These experimental efforts have lead so far to 
the following picture^: for strongly disordered packings (for example by considering mixtures of grains of very different sizes, 
or irregular grains such as natural sand), the response profile on large length scales shows a single broad peak [69,116]. This 
single hump response function was also observed in numerical simulation of 2D polygonal grains packing [101]. 

For well ordered packings however, the two peaks structure is rather convincingly observed in two dimensions (a 'ring' or 
three peaks in 3D) [69,103,15]. These hyperbolic features are also in agreement with a recent numerical simulation on (rather 
small) isostatic assemblies of frictionless grains [138,78], with the work reported in [63,28] where beads are arranged on a regular 
triangular lattice but where disorder is introduced by a finite friction coefficient. Similarly, two-peak response can be observed 
for strongly anisotropic elastic networks [73,110]. 

Obtaining a precise experimental determination of the linear response function that can be quantitatively compared to 
theoretical models is rather difficult: the perturbation must be small enough not to disrupt the packing, but also large enough 
to lead to a measurable signal. A very sensitive technique, based on a lock-in detection of an oscillating perturbation, has 
allowed one to obtain precise and reproducible results [116,124]. One finds unambiguously that the response function G(r,H) 
to a perturbation located at ro = zq = is in the case of sand layers single-peaked, with a width growing as the height H of the 
layer. More precisely, the response function for different heights can be rescaled onto a unique curve by plotting H^G{r, H) as 
a function of r/H. The factor is expected from force conservation: the integral of the response function over the bottom 
plate must be equal to the overload force F for the total force to balance. One also finds that, like for sand-piles, the pressure 
response profile depends on the way the granular layer was prepared its 'history': the value of the maximum of this response is 
roughly twice smaller (and its width twice larger) for a dense packing than for a loose one (see Fig. 6). The quantitative shape 



^Note however that a purely 'diffusive' response function, scaling as VH as predicted by the 'scalar' model, was reported in 
[126] for a very special 'brick' packing. 
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of the experimental response function cannot be accounted for by a simple isotropic elasticity theory (sec also the discussion 
in [147]), but in the present geometry, anisotropic elasticity leaves three extra adimensional constants that can give a large 
freedom to fit the data (see [110]). 

Therefore, the most precise experiments on strongly disordered packings seems to favor a single peak, elastic like response, 
rather than a double peak response that only appears in special conditions (ordered packings, or frictionless grains). We 
therefore need to understand in details the role of randomness in hyperbolic wave equations, or in more physical terms, the 
large scale consequence of the fact that force chains are 'scattered' by packing irregulaxities. Is this sufficient to convert a two 
peak response function into a single peak, elastic like function? This is the subject of the following sections. 



VII. FORCE CHAINS SCATTERING I: WEAK DISORDER LIMIT 



A. A stochastic wave equation 

Provided that local conservation laws (those arising from mechanical equilibrium) arc obeyed, many local rules for force 
transmission are a priori compatible with the existence of contacts among rigid particles [63,131]. Therefore, even if there is a 
definite mean relationship among stresses at the meso-scale, one can expect randomness in the local transmission coefficients. 
The simplest model for this and other sources of randomness is to introduce a randomly varying 'speed of light' Co- This could 
describe the fact that, for example, the parameter p in the model of figure 5 can vary from grain to grain. In this situation the 
two rays arc still symmetric around the vertical direction, but with a random opening angle. The situation where the bisecting 
line itself is random (i.e. when both the above opening angles tp± are fluctuating) will be alluded to below. 

This suggests the following stochastic wave equation for stress propagation in two dimensions: 

dttCTtt = dxx [co(l + v{x,t)) crtt] (71) 

where the random noise v is assumed to be correlated as {v{x,t)v{x' ,t')) = a^g^ix — x')gt{t — t'). The correlation lengths 
Ix and It are again kept finite, and of the same order of magnitude. In Fourier transform, this relation can also be written 
{v{k,t)v{k' , t')) = 2na^5{k+k')gx{k)gt{t—t'). It turns out that the final shape of the averaged response function depends on the 
sign of §x(A). In section III we implicitly made the choice gxik) = 1, which corresponds to: gx{x = 0) = 1/ a and gx{x > 0) = 0. 
We will keep this choice for the following calculations, but note that another form for gx could lead to gx{A)) < 0. 
In the following, att will be again denoted by w. After a Fourier transform along x-axis, we get, from equation (71) 

{dtt + ($k'^)w = dtp - cofc^ J ^^(S; - 1, t) (72) 

Note that the 'source' term of this equation is now dtp rather than p itself. Therefore, if we call G the Green function (or 
propagator) of this equation G = {5w/5dtp}; the response function R = (Sw/Sp) of our system is now actually the time 
derivative of G: R{k,t) = dtG{k,t). 
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The noiseless propagator Go is the solution of the ordinary wave equation {dtt + Cok^)Go{k,t — t') = d{t — t') and can be 
easily calculated: 

Go{k, t) = [e''^"** - 6"*'^"'=*] d{t) (73) 

which leads to the response function Rq 

Ro{x, t) = ^[6{x- cot) + 6{x + cot)] e{t) (74) 

This last equation sums up one of the major results of the hyperbolic approach of [18,145,146]: in two dimensions, stress 

propagates along two characteristic rays. (Note that the corresponding response function in three dimensions reads i?o(a;,t) oc 
(cgi^ — x^)^^^^ for |a::| < Cot and zero otherwise [18]). A relevant question is now to ask how these rays survive in the presence 
of disorder. We will show that for weak disorder, the 5-peaks acquire a finite (diffusive) width, and that the 'speed of light' is 
renormalized to a lower value. Not surprisingly, the effect of disorder can be described by an 'optical index' n > 1. For a strong 
disorder, however, we find (within a Gaussian approximation for the noise v) that the speed of light vanishes and becomes 
imaginary. The 'propagative' nature of the stress transmission disappears and the system might behave more like an elastic 
body, in a sense clarified below. 



B. Calculation of the averaged response function 

One can again use Novikov's theorem in the present case if the noise is Gaussian and short range correlated in time. However, 
the same results are again obtained within the diagrammatic approach explained in section III, which can be easily transposed 

to the present case, and is more general. The propagator G is a now represented as a line, the source dtp a cross and the vertex 
meaning — Cq/c^ J Within the MCA, the self-consistent equation (analogous to equation (38) in the scalar case) is: 



fdt' 
Jo 



{dtt + Cok^)H{k,t) = 6{t)+ / dt'T,uoA(k,t')H{k,t-t') (75) 



where H is defined by G{k,t) — H{k,t)8{t), and the self-energy ^MCA given BS 

^uoA{k,t-t') = 2nctc7^k^ j ^q'gt{t-t')g4k-q)H{q,t-t') (76) 

Equation (75) can be solved using a standard Laplace transform along the f-axis {E is the Laplace variable). Using the 
fact that H{k,T) = t in the limit where t — > 0, we find, for small k,E (corresponding to scales L such that ix, £t <C L): 
H-\k, E) = E'^+I3E + c\k'^, where 

4(.)=cg-^(l-f)+0(.^) (77) 
/3(fc) = £o^AAA+o(fc^) (78) 

We notice here that in the limit It — > 0, the effect of the randomness completely disappears, as in the scalar model with the 

Ito convention. (Technically, this is duo to the fact that G{k,t = 0) = in the present problem). In order to calculate the 
inverse Laplace transform, we need to know the roots of the equation H~^{k, E) = 0. This leads to several phases, depending 
on the strength of the disorder. 

• The weak disorder limit. 

For a weak disorder, c%{k) is always positive. Wo can then define c_r — CR{k = 0). As we will show now, cr, is the shifted 
'cone' angle along which stress propagates asymptotically, cr is a decreasing function of a, meaning that the peaks of the 
response function get closer together as the disorder increases ^ . For a critical value a = ctc, cr vanishes, and becomes imaginary 
for stronger disorder. 



^As a technical remark, let us note that if gt = gx, the problem is symmetric in the change x —* t, c^{x,t) l/(:^{x,t). It 
thus looks as if the cone should both narrow or widen, depending on the arbitrary choice of x and t. There is however no 
contradiction with the above calculation since we assumed that v has zero mean, while 1/(1 + v) — 1 has a positive mean value, 
of order 
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In the limit of large t, the propagator reads: 



G{k, t) = ^ sin [cRkt{l + a\k\)] e-^'=^*6l(t) 

CrK 



where the following constants have been introduced 



7 



I3{k) _ a^A^i^ 



2fc2 18 

From equation (79), the response function 7?, in the limit of small k and large t, is given by: 

R{k,t) = cos [cRkt{l + a\k\)] e"^'"'*6'(t) 

or in the real space, 



R{x,t) 



where the scaling variables measuring distances relative to the two peaks, are defined by 



(79) 

(80) 
(81) 

(82) 
(83) 




FIG. 7. Response function for weak disorder {a/ac ^ 0.32). The two curves have been rescaled by the factor 2 [47r|7|t]^''^. 
The main graph shows the general doublc-pcakod shape of the response of the system when subjected to a peaked overload at 
X = 0, t = 0. The inset gives details the right-hand peak as a function of the scaling variable Note the asymmetry (for 
gx(A) > 0), compatible with the results found in [63]. Note also that the curve becomes negative around ^- = 2. 



X ± CRt 

7^ 



(84) 



and where 7 = 'y—iCRa, b = e'"^'''. $ is the standard error function. Figure 7 shows R as given by expression (83). Interestingly, 

this propagator not only has a finite diffusive width oc ^/i, but is also asymmetric around its maxima. Surprisingly, and in sharp 
contrast to the scalar case discussed above, the response function becomes negative in certain intervals (although its integral 



*Note that the sign of a is dictated by the sign of 5x(A) 
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is of course equal to one because of weight conservation). This means that pushing on a given point actually reduces the 
downward pressure on certain points. This can be interpreted as some kind of arching; increasing the shear stress does affect 
the propagation of the vertical stress, and may indeed lead to a reduction in its local value which is redistributed elsewhere. 
The un-averaged response function therefore necessarily takes negative (and in fact rather large [41]) values. This is actually 
of crucial importance since this is the source of fundamental 'fragility' of granular matter to external perturbations. Suppose 
indeed that as a result of the perturbation, a grain receives a negative force larger than the preexisting vertical pressure. This 
grain will then move and a local rearrangement of contacts will occur, inducing a variation of co{x,t) as to reduce the cause 
of the instability. Thus, the stochastic wave equation implicitly demands rules similax to those introduced, e.g. in [40]. The 
present model, which is purely static, does not say what to do when a local rearrangement occurs, but certainly suggests that 
small perturbations should induce such rearrangements. 

It is interesting to note that this response function was numerically measured in [63] ; its shape is compatible with the above 
expression; in particular, the two peaks were found to be asymmetric with a longer 'tail' extending inwards, as we obtain here. 
Note however that for gx{J^) < 0, the shape of the peaks is reversed: the small dips are located inside the peaks and the longer 
tail extends outwards. In a more recent work [28], both the diffusive spreading and the renormalisation of the wave velocity 
have been observed. 

• Critical disorder: The wave/diffusion transition. 

When the disorder is so strong that cr just vanishes, the roots of H~^{k,E) = change nature, and so does the response 
function R. The two peaks of the previous expression for R merge together, while the width becomes anomalously large 
(oc t^^^). In the asymptotic, large t, regime we obtain: 

R{k, t) = e(t) cos [X\kf''^t\ e"^*'* (85) 

where the new constant A is defined by A = co •\/3/2A and 7 = Co^t /3. 

• Strong disorder: The pseudo-elastic regime. 

For larger disorder still, one finds, within the MCA (which is exact for a Gaussian, uncorrelated noise), that the renormalized 
value of Co, c%,. becomes negative. Upon a rescaling of a; as x = x/{icR), the effective equation on {5att) then becomes, on large 
length scales, Poisson's equation: 

{6(Jtt) = dt {5p) (86) 

which means that the stress propagation becomes somewhat similar to that in an elastic body, where stresses obey an elliptic 
equation of similar type [89]. In particular, the cone structure of stress propagation, whicii is associated with the underlying, 
hyperbolic, wave equation finally disappears; the average response to a localized perturbation becomes a broad 'bump' of width 
comparable to the height of the pile. However, the above transition is possibly an artifact, due to the fact that v is taken to be 
Gaussian, which strictly speaking is not allowed, since the local value of should always be positive. One can show for some 
other problems of the same type that a similar transition is artificially induced by the Gaussian approximation when it cannot 
really exist on physical grounds. We need to address the strong disorder limit with different tools in order to discuss the large 
scale nature of the response function in this case: this will be discussed below. 



C. Generalized wave equations 

It is tempting to generalize equations (58, 59) and write the most general linear equations governing the propagation of the 
forces which are compatible with the (local) conservation rules. These equations were first written by de Gennes [70]: 

dtFt + dx [n'{x, t)Fx + v'{x, t)Ft] = p (87) 
dtFx + dx[ri{x,t)Ft + v{x,t)Fx] =0 (88) 

Note that the terms v, v' break the symmetry x —x, and exist whenever the transfer rules in the three leg model are not 
symmetrical. Equations (87, 88) describe a situation where not only the opening angle of the cone can vary in space, but also 
its average orientation. 

The same analytical techniques as above can still be used. We shall only discuss some special cases ^: 



^To lowest order in perturbation theory, the case where disorder is present in the four terms rj,'!]' ,v, v' simultaneously is very 
simply obtained by adding the corrections induced by each term taken individually. 
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o Let us first set v — v' = and consider the case where r;' is random, and r; fixed (and equal to Cq). Taking ri'{x,t) = 
rjo (1 + v{x,t)) with the noise v as above, one finds that the rcnormalized value of rj' is: 

= »7o 1 7, (89) 



Now, on large length scales, one must recover the continuum equilibrium equations for the stress tensor, equations (60, 61). 
The condition of zero torque requires that the stress tensor is symmetric, and thus one must set 77^ = 1, which imposes a 
relation between 770 and the amplitude of the noise a. Note that beyond a certain value of a, this relation can no longer be 
satisfied with a real 770 . This again means that the packing is unstable mechanically and will rearrange so as to reduce the 
disorder. 

o Another interesting class of models, which one can call '/u models', is such that: 77 = 00,77' = 1, but n{x,t) = cov{x,t) and 
/Lt' = or vice-versa. These two cases yield identical results, namely, in the large t limit: 

R{k, t) = cos (cofct) e"^*''6i(t) (90) 
i?4fc,t) = -ICO sin (cofct) e"^*^*6»(t) (91) 

where 7 = — . Note that in these cases, the response peaks acquire a finite diffusive width oc \/t, but the angle of the 
information cone is unaffected by the disorder (i.e. co is not renormalized) . 

o Finally, there are special 'symmetry' conditions where the equations can be decoupled and reduced to two 'scalar' models. 
We will refer to this case as the 'double scalar' model. This occurs when n = fj,' = coVi(x,t) and 77' = 77/co = 1 + V2{x,t) where 
vi,V2 are two different sets of noise. Let us define cr± = co^t ± Fx, x± = x^ cot and v± = vi i: V2, we then obtain: 

dtCT+ = Cop - Co dx^ [v+a+] (92) 
dta- = Cop - Co dx_ [v-o-] (93) 

showing that a+ and a- decouple, each propagating along two rays, of 'velocity' ±co, plus a small noise v± which, as in the 
scalar case, generates a nonzero diffusion constant. The response functions for crtt and a^t are thus again made of two diffusive 
peaks of width oc \/t, centered in x = ±cot. Note that by construction, this special form of disorder does not lead to negative 
vertical stresses. 

A physically relevant question is to know how local stresses are distributed. We have seen above that within a scalar 
approach, an exponential-like distribution (possibly of the type exp—w^, with /3 > 1) is expected [93,46]. One can wonder 
whether this exponential distribution survives within a tensorial description, and what happens for very small stresses w — > 0. 
Unfortunately, the full distribution can only be computed analytically for the 'double scalar' model; but numerical results have 
also been obtained for the random BCC model [41] and confirm that the exponential tail holds in the strong disorder limit. 

In the 'double scalar' limit, the histogram of the stress distribution is obtained trivially by noting that since cr+ = w\ and 
cr_ = W2 travel along different paths, they are independent random variables. Taking co to be unity for simplicity, one thus 
finds: 

P{au) = J d^ij dw2P*{wi)P*{w2)5{au - !^^1±^) (94) 
PicTxt) = JdmJ dw2P*{wi)P*{w2)S{a^t - ^1 - ) (95) 

where P* is the distribution of weight pertaining to the scalar case, which, as mentioned above, depends on the specific form 
of the local disorder and on the discretisation procedure. In the strong disorder case which leads to equation (8) [in the case 
N = 2], we thus find that P{crtt) is still decaying exponentially (it is actually a F distribution of parameter 2A'^), although its 
variance is reduced by a factor 2. For N = 2, one simply gets 

P{<ytt) = ^af^e-^"^" (96) 
P(cr.,)=(|a.,| + i)e-^l'^-l (97) 

The interest of the 'double scalar' limit of the hyperbolic model is to show that the exponential tail Eq. (8) has indeed a certain 
degree of universality, and is not restricted to the scalar model. 
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middle of the top surface. To get this picture birefringent grains between inverse circular polarizers have been used, and the 
intensity difference after and before the overloading has been computed. This picture was obtained by R.P. Behringer and J. 
Gang. 



VIII. FORCE CHAINS SCATTERING II: STRONG DISORDER LIMIT 
A. Introduction and numerical results 

Prom pictures of photo-elastic grains, the network of inter-particle forces propagating as a responses to a localized pressure 

was extracted [43,69] (see Figure 8). An interpretation of such a picture can be given in terms of linear force chains which tend 
to split upon meeting vacancies or packing defects [24]. 

As reviewed in the previous section, one can investigate perturbatively the role of disorder on hyperbolic equations. In this 
case, the two peak structure of the response function is preserved on large length scales, although the peaJss are diffusively 
broadened. This result is in qualitative agreement with the numerical simulations [63,78,28]. An uncontrolled extrapolation to 
strong disorder however suggests that the large scale equations might become elliptic. 

In order to investigate more quantitatively the strongly disordered regime where force chains split and merge, one can study 
[24] the following model. If one of the force chains meets a defect (randomly distributed in space), we split it into two new ones 
at random angle, which then propagate until another defect (or the boundary) is reached. More precisely, a chain carrying a 
force / in the direction ft splits when meeting a defect into two forces /i, /2 in the directions ni, n2 - 'A process'. The two 
angles 6i and 62 (between ft and ni and n2 respectively) are uniformly chosen between and some maxinnim splitting angle 
9m- The local mechanical equilibrium imposes that the intensities /i and /2 arc such that fft — fifti + f2n2- Sometimes, two 
(or more) force chains meet at the same defect - 'Y process'. In this case, we make them merge together. It is important to 
note that the positions of the defects are fixed before starting the computation of the forces. This idea of a frozen disorder is 
consistent with the experimental observation that when the local overload is added on the top of the system, the forces are 
transmitted along the chains originally created during the building of the packing. In other words, as long as the applied force 
is not too large and compatible with the pre-existing network of force chains, the geometry of the packing, and in particular 
the contacts between grains, remains the same (but see the discussion in section IX). 

With these rules, realistic force networks can be created numerically [24]. After averaging over many statistically identical 
samples, one can obtain stress profiles for different heights H. One finds that, as H increases, the vertical pressure response 
profile evolves continuously from two well defined peaks to a single broad one. It means that the hyperbolic behaviour is 
progressively erased by multiple scattering. The width of the single peak is found to scale like H; the scaling function is 
furthermore surprisingly close to the pure elastic response of a semi-infinite two-dimensional medium. However, the numerical 
simulation was restricted to rather shallow (small H) samples. 
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B. A Boltzmann description of force chain splitting 



In order to understand analytically the above numerical results, we write a Boltzmann equation for the probability density 
P{f,n,r) of finding an oriented force chain of intensity / in the direction n around the point r [24]. A very important point 
here is that force chains can be oriented in reference to the boundary conditions (see the discussion in [34,138,24]). 

For simplicity, we first neglect the chain 'merging' process which leads to a more complicated non linear Boltzmann equation 
(its influence is not fully understood yet and will be discussed below). We also assume that the splitting is symmetric, i.e that 
ft- Hi = n • n2 = cos 6, so that /i = /2 = //2 cos 6. Assuming a uniform density of defects, the probability distribution P{f, ft, r) 
obeys the following general equation: 

P{f,n,f+ndr) = (1 - ^)P{f,fi,r) +2^ Jdn' Jdf P{f' ,n' ,f)^{n' ,n) 5(^f - (98) 

where A is equal to the 'mean free path' of force chains, and is of order l/(pda'*~^) in dimension d. The above equation means 
the following: since a chain of grains can only transport a force parallel to itself [35], the direction of the force n also gives the 
local direction of the chain. Between r and r + ftdr, the chain can either carry on undisturbed, or be scattered. The second 
term on the right hand side therefore gives the probability that a force chain initially in direction n' is scattered in direction n. 
This occurs with a probability ^^{n', n), where is the scattering cross section, which we will assume to depend only on the 
scattering angle 0, for example a uniform distribution between —9m and +9m- The J-function ensures force conservation and 
the factor two comes from the counting of the two possible outgoing force chains. Let us now multiply equation (98) by / and 
integrate over /. This leads to an equation for the local average force per unit volume in the direction n, that we will denote 
F(n,r). This equation reads: 



\n-VrF{n,r) = -F(n,r) + dn I ' ^(n ,n) + -n-Fo(r), 

J n ■ n a 



(99) 



where we have added the contribution of an external body force density Fo{r), and a is the size of a defect or of a grain. 
This equation is the so-called Scliwarscliild-Milne equation for radiative transfer, describing the evolution of light intensity in 
a turbid medium [142]. We now introduce some angular averages of F{n,r) that have an immediate physical interpretation: 

p{r) = ajdQ.F{n,r) (100) 
J^{r) = ajdQ naF{n,r) (101) 
(^a/sir) = ad dQ nctUg F{n,r), (102) 



where JdCl is the normalized integral over the unit sphere, that is introduced for a correct interpretation of cr (see equation 
(106) below). As will be clear from the following, J is the local average force chain intensity per unit surface, a is the stress 
tensor. Since = 1, one finds that Tr a = dp, and therefore p is the isostatic pressure. Note that J is not the average local 
force, which is always zero in equilibrium. The fact that J ^ comes from the possibility of orienting the force chains. 

We now integrate over the unit sphere equation (99) after multiplying it by different powers of ria- Using the fact that 
'^{fi' , n) only depends on n • n', a direct integration leads to: 

XV ■ J = (ao - l)p, (103) 

where ao is called the 'albedo' in the context of light scattering [142], and reads: 



ao= fdn '^^f'^ > 1. (104) 
J n-n' 

A second set of equations can be obtained by multiplying by ria and integrating. Using the fact that jdn'^{n\n) = 1, it is 
easy to show that 

Idnn \ I ' = n . (105) 
J n ■ n' 

Therefore, the resulting equation is nothing but the usual mechanical equilibrium relation: 

Vaa„,a = Foc. (106) 

This relation in fact only reflects the local balance of forces chains. Now we multiply equation (99) by nanp and again integrate. 
A priori, this introduces a new three index tensor. In order to close the equation, we now make the an assumption that is 
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usually made in the context of light diffusion, that on large scales the force intensity becomes nearly isotropic [142]. In this 
case, it is justified to expand F{n, f) in angular harmonics and to keep only the first terms: 

aF{n,T)=p{f)+dn- J{r) + ... (107) 

Within this expansion, we finally obtain a 'constitutive' relation between aafj and the vector J. We find: 



< dao — an 



1) 



(Vc J/3 + V/3 Ja) , (108) 



and 



Ad 



f do2— gp _ -I \ 
d-1 



2dKl Vc Ja + ( dJCi - "° , 1 V- J 

(d- l)(oo - 1) 



(109) 



with K\ = l/d(2 + d) and 02 = fdn{n • n')^(n',n). Rather surprisingly, these equations have exactly the canonical form of 
elasticity theory, provided one identifies the vector J with the local displacement, up to a rrmltiplicativc factor. 

The above stress equations are rather non-trivial because no displacement field is introduced in the above derivation (nor in 
the numerical model): elastic- like equations are found in a stress-only model. The basic assumption is the existence of local 
force chains, which have a well defined identity over several grain sizes a, such that a ^ A: this is the condition under which 
the above Boltzmann description of the force chain scattering is justified. 

Since the above equations are formally identical to those of classical elasticity, one can show that V^p = 0, and V'j — 
[89]. One can therefore in principle compute the response function G(f) to a localized force at fo = in the z direction, which 
is given by the standard (one peak) elastic Green's function. But note however that although the above equations are formally 
those of classical elasticity, there is one crucial difference coming from the fact that (cra/3) and J are not independent since they 
are both related to the same underlying quantity F{n,r). This is a very important difTerence, which appears, for example in 
the choice of boundary conditions on B that determines P{f, n, 



C. The role of chain merging 

We have mentioned above the fact that the numerical simulation of the full chain scattering model (including both splitting 

and merging) was restricted to small depths. Similarly, neglecting chain merging altogether cannot be correct on large length 
scales, since the number of force chains would diverge, leading to an infinite number of force chains with infinitesimal intensity. 
As shown in [132], the above linearized theory is in fact unstable, and chain merging play a crucial role to make the theory 
mathematically consistent. 

On the other hand, chain merging leads to a non-linear Boltzmann equation which in the general case cannot be solved. 
Only special cases have, up to now, been amenable to an analytic treatment. For example, if one insists that the force chains 
can only take six directions a 60° degrees, with 120° degree chain splitting and chain merging, the amplitude of the force chains 
is constant, simply because 2 cos 60 =1. In this case, the full probability distribution P{f,n,r) boils down to six functions 
Pn{r), describing the probability for a force chain to be in one of the six available directions. The result of the full analysis is 
that, quite surprisingly, the response function evolves back to a two-peak, hyperbolic structure on large length scales! Whether 
this is duo to the particular structure of the model is not yet settled; preliminary results on the eight-fold model [26] suggests 
that there might actually be a transition between a hyperbolic response for sufficiently anisotropic scattering and an elliptic 
response for isotropic situations, (see [133] for a nice recent discussion). The existence of such a transition would perhaps, as 
we discuss now, allow one to reconcile the apparently contradictory experimental, numerical and theoretical results. 



IX. STATICS OF GRANULAR MATERIALS: CONCLUDING REMARKS AND OPEN QUESTIONS 

Let us try to summarize the theoretical situation as follows. If the grains are frictiorilcss, then packings are generically 
isostatic. If these packings are furthermore 'sufficiently' anisotropic, then the construction of the stress everywhere in the 
packing can be performed by 'propagating' the stress from one boundary towards the interior of the packing, as one expects for 
hyperbolic equations. This situation is well captured by the three-leg model introduced above, that indeed leads to hyperbolic 
equations on large length scales. Correspondingly, numerical simulations of the response function in anisotropic packings of 
frictionless beads indeed display the expected diffusively broadened two-peak structure, at least for the modest sizes that were 
simulated [78]. 

However, not all isostatic packings are characterized by a two-peak response function. For example, using isotropic isostatic 
packings built by J.N. Roux, Ph. Claudin and A. Ayadim [1] have found that the average response function is a single hump, 
qualitatively similar to the elliptic-like response observed experimentally and to the result of the chain splitting model. An 
example of such packings is shown in Fig. 9, together with the local force chains. This picture makes it obvious that the 
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FIG. 9. Force chains in an isotropic, isostatic packing of frictionless disks. The average response function in this case has a 
single peak. Courtesy of Ph. Claudin. 



force chains are strongly scattered, and lose their 'coherence', which makes it indeed plausible that the two-peak structure 

is destroyed. Again, some degree of anisotropy seems to be required to preserve hyperbolicity. Note that in the numerical 
determination of the response function, the applied force perturbation 5f is infinitesimal and does not induce any rearrangement 
of the initial structure. 

This remark is important since i/ rearrangements are allowed, the response has been argued in [118] to be, on general grounds, 
elliptic. Moreover, one expects that for any small perturbation, a sufficiently large assembly of frictionless grains will rearrange 
[44]. Therefore, the order in which the limits H ^ oo and 5/^0 arc taken is physically relevant [119,79], and one might 
expect that even anisotropic structures such as those studied in [78] should destabilize under a small perturbation for large 
enough H. It is not clear which of the limits is relevant in the experimental conditions of [116]; even if the applied force is very 
small and much care has been devoted to perturb as weakly as possible the packing, it might well be that these experiments 
are not in the infinitesimal limit. 

Conversely, is isostaticity needed to obtain hyperbolicity? The numerical simulations of [63,28], where an hyperbolic response 
function for anisotropic ordered packing of grains with friction is observed, show that this is not the case. Similarly, the 
experimental determination of hyperbolic like response function in ordered lattices [103,69] shows that, as suggested by the 
above theoretical analysis, weak disorder does not suppress the hyperbohc nature of the force propagation. It would be 
extremely important, in this context, to confirm theoretically the existence of a true 'hyperbolic-elliptic' phase transition in 
simplified models of force chain scattering. 

More formal approaches have also been advocated in [62,2] to try to establish some closure relations between the components 
of the (coarse-grained) stress tensor starting from the equilibrium force balance at the grain level. Tkachenko and Witten [138] 
have insisted on the equivalence between any macroscopic linear relation between the components of the stress tensor and the 
existence of 'fioppy' modes, i.e. zero energy deformation modes (see also [102]). Such floppy modes can exist in certain elastic 
spring networks, such as a square lattice which can deform at no cost along its diagonal. One can indeed show that the large 
scale elasticity equations in this case are characterized by coefficients such that the marginal condition = t (see Eq. (54) 
above) is satisfied, so that the response function becomes the sum of two delta peaks [110]. Following this line of thought, 
a natural conjecture is that if a system has some extended floppy modes, then its large scale response will be hyperbolic. 
Conversely, if floppy modes are all localized, then the response is locally hyperbolic, but the strong force chain scattering 
disrupts the long range propagation and the response becomes elliptic like, as suggested by the force chain scattering model 
discussed above. We hope that these issues will be clarified in the near future; a particularly important point is to establish 
precisely, starting from a microscopic description a la Edwards, which kind of large scale static equation emerge and under 
what conditions it is hyperbolic/elliptic, and the value of the parameters entering these equations. 
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X. GLASSY DYNAMICS IN GRANULAR MEDIA: A BRIEF SURVEY 



A. Slow compaction 

We have spent quite a long time on the static properties of granular media. The dynamics of weakly tapped assemblies 
of grains is also an extremely interesting subject, in particular in relation with the dynamics of other glassy systems, such 
as supcr-coolod liquids, colloidal glasses and foams. The most obvious experiment (with important industrial applications) 
is that of compaction under tapping. The basic experiment consists in studying the volume occupied by a large number of 
grains when the container is 'tapped', i.e. subject to a periodic acceleration of a certain amplitude. The packing, initially very 
loose, progressively compacts. However, the compaction process slows down with time, and the decay of the volume is very 
far from being a simple exponential. Experiments have shown that the decay of the volume towards its asymptotic value can 
be satisfactorily fitted by an inverse logarithm of time, or by a stretched exponential, both forms being typical of relaxation 
processes in glasses, where one can also study the time dependence of the volume (or of the energy) as the sample relaxes after 
a temperature quench. The parameters of these fits depend on the tapping amplitude - stronger tapping obviously leads to a 
faster compaction. 

More complex experimental protocols have also been tested. For example, one can change, in the course of the experiment, 
the amplitude of the tapping and reveal interesting memory effects, again similar to those found in glasses and spin-glasses. 
The now clEissic experiment [108] is to start from a loose sample and incrcEise slowly the tapping amplitude F, in such a way 
that for each tapping amplitude the density p appears to reach a saturation value. One finds that p(F) increases with F. When 
F is reduced back to zero, the density keeps a high value, revealing a kind of irreversibility that also appears in spin-glasses 
under a magnetic field: when the temperature is increased the (zero field cooled) magnetization increases, but does not follow 
the same path on the way back. The temperature at which the two branches separate is the spin-glass transition temperature 
(which, if defined in this way, weakly depends on the heating/cooling rate). Similarly, there appears to be a tapping amplitude 
beyond which the two density branches meet. The high magnetization (field cooled) branch, as the high density branch in the 
granular system, is reversible. As in spin-glasses, the low density branch is in fact out-of-equilibrium, but the convergence of 
the density (magnetization) towards its equilibrium value is much too slow to be observed. 

In the first stage of another type of experiment [83], aimed at revealing 'memory effects', one taps the system with three 
different amplitudes - say weak, moderate and strong - during a time chosen such as to reach a certain density, identical in the 
three cases. In the second stage of the experiment, the tapping amplitude is then chosen to be moderate, and the density just 
after the amplitude 'jump' is observed. If the state of the system was only described by the density, the evolution of the density 
after the jump should be identical for all three situations, and follow the 'moderate' curve, which is taken as the reference. This 
is not the case: the weakly tapped system first has to dilate before it is able to resume its compaction, whereas the strongly 
tapped system compacts faster than the reference system just after the jump. This shows that the configurations reached 
under stronger tapping are in a sense easier to compact further. This experiment therefore indicates that some further 'hidden' 
observables arc needed to describe the large scale evolution of the system. We shall expand on this below, but want to note 
that a very similar effect is known in glasses as the Kovacs (or memory) effect. In this case, one prepares the system at a given 
temperature T2 < Ti and waits until the volume has reached the equilibrium volume at Ti. Then one raises the temperature 
to Ti. If the volume was the only relevant quantity, one should not see any further evolution since this volume is already at its 
equilibrium value. Again, this is not what Kovacs first observed in polymeric glasses [84]: the volume has to increase back to 
a maximum before decreasing again towards its equilibrium value. A similar effect was also reported in numerical simulations 
of spin-glasses [9] and of Lennard- Jones systems [12], and several simple models are now known to reproduce this effect (see 
below). 

We also want to mention the very interesting experiment of d'Anna et al. [48], where a torsion oscillator is immersed in a 
vibrated granular assembly. The observable is the angle 9 made by the oscillator with an arbitrary axis. The results of [48] are 
that (a) the angle 9 performs a random walk in time: {[9{t + T) — 6{t)]^) = D{r)T and (b) the angular diffusion constant D{r) 
appears to vanish as the amplitude of the tapping F goes to zero in a way that recalls the 'super-activated' Vogel-Fulcher law 
in glasses. 

B. Self-inhibitory dynamics and dynamical heterogeneities 

1. Non exponential relaxation 

We shall call p* the maximum compacity that can be reached in a tapping experiment. For F — + 0, this corresponds to 
the so-called 'random close packing' configuration. [It is not entirely clear how meaningful this notion really is, but from an 
empirical point of view, any reasonable extrapolation of the accessible dynamics seems to converge to a density which is not 
the HCP density (in the case of identical spheres).] Correspondingly, we will call the difference between p* and the average 
density p the 'free volume' fraction, The simplest relaxation equation for $ is obviously a rate equation: 
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- = (110) 

If the decay rate 7 is independent of the free-volume itself, the decay of $ is obviously a single exponential. However, it is 
intuitively clear that the dynamics is self inhibitory, in the sense that it is the free volume itself that allows further compaction. 
Therefore, one expects that 7 vanishes as $ — > 0. Assuming a power-law dependence 7 ~ with /3 > 0, one obtains a 
power-law relaxation for long times: 

^^(t + to)-'^''. (Ill) 

Let us assume a simple model where particles have a volume u, and mobile holes have fixed 'quantum' volume vo- For a particle 
to be able to move, we require that n* holes meet in a volume v, such that the volume of one particle is liberated, i.e. n*vo = v. 
If the dynamics is sufficiently chaotic, it is reasonable to assume a Poisson distribution for the holes. Therefore, the probability 
for a particle to move is simply , leading to /3 = n* . If the number of holes n* needed to move one particle is large, one 
can approximate the long time behaviour of 3> as a logarithm: 

, (112) 

a form often advocated to fit the slow relaxation of the volume in glasses or granular media [135,108] (but see also [111], where 
the role of convection is discussed, and [137] for a simple soluble model). 



2. Dynamical heterogeneities 

If one assumes that holes cannot spontaneously appear, it is clear that the dynamics is necessarily spatially heterogeneous: 
regions rich in holes, whore dynamics is locally fast, appear only to the detriment of regions poor in holes, whore the system 
is jammed. This argument can be somewhat formalized to suggest that the geometry of the 'fast' objects is non trivial. Call 
3> + <i>{r,t) the coarse-grained density of free volume ('holes') around point fat time t, such that the space average of <f){r,t) 
is equal to 0. [Note that all 'voids' do not necessarily contribute to the free volume). Far from the boundaries of the sample 
where the holes can disappear, one can write a Langevin equation for (j>{f,t) which reads [51]: 

|^=DVV + V-[V* + 0eW], (113) 

where D is the (fast) diffusion constant of the holes, ^ is a Gaussian white noise in space and time of variance equal to D, and 
we have neglected the interaction between the holes. In the above equation, we have supposed that the free-volume is locally 
conserved. Assuming that (j) is small leads to the result that itself is a Gaussian field with a correlation function given (in 
three dimensions) by:" 

(Hr,tmrM + r)) = | (^^^exp (-^^) (r > 0). (114) 

Note that in the limit r ^ 0, the field <j) is dclta-corrolatcd in space; some short scale cut-off of the order of the grain size a 
would be needed to regularize this behaviour. Now, if one insists that the material particles (or grains) can only move if the 
local density of holes is larger than a certain threshold $c, the diffusion equation for, say, the density p{f,t) of slow tracer 
particles reads (see also [94] for similar ideas): 

^=DoV- [e($ + ct>- *c)Vp] . (115) 

Although a more careful solution of the coupled equations (113, 115) is required to confirm the following conclusions, a qualitative 
analysis of the problem leads to the following picture for glassy dynamics, which relates slow, self-inhibitory dynamics and 
dynamical heterogeneities: 

• The probability that an elementary region of space is active is, for small $ <^ $c, proportional to exp(— C$c/^)) where 
C is a numerical coefficient. Therefore, the large scale diffusion coefficient Dr of the particles vanishes as: 

Dr^ Doexp(^-^^ =Doexp(^-^^^ , (116) 

i.e. a la Vogel-Fulcher. This argument is actually the analogue of the classic free-volume argument leading to the 
Vogel-Fulcher law in glasses, and extended to granular media by Boutreux and De Gennes [27]. 



If one considers that 'holes' can be converted into voids and vice-versa, such that 4> is only conserved on average, the statistics 
of the 4> field is altered and becomes that of an elastic membrane subject to thermal fluctuations. 
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• The fast regions arc delimited by the contour lines of a correlated Gaussian field. This has interesting consequences: for 
example, one expects, in the short-ranged conserved case, the active regions to be lattice animals of fractal dimension 
d/ = 2 in three dimension (and dj ~ 1.56 in two dimensions). The value df = 2 turns out to be in agreement with the 
experimental determination of Weeks et al. on dense, three dimensional colloidal glasses [144]. It would be interesting 
to measure other geometrical characteristics of the fast clusters to confirm that these are indeed lattice animals (see e.g. 
the discussion in [88]). Interestingly, the connection between glassy dynamics and the contour lines of a dynamic random 
field was also suggested in [36], using rather difi'erent arguments. 

It would be most interesting to make these statements more precise in the context of a specific model. Very promising steps 
in that directions have been made recently in dynamically constrained models [67,139,117]. 

As emphasized by Boutreux and De Gennes, a Vogel-Fulcher law such as Eq. (116) also leads to a logarithmic relaxation of 
the density. Indeed, writing: 

( C^l\ ^ 

= -7oexp ^ (117) 



dt ' " V * 

leads at large times to: 



(118) 

log t 



3. Another point of view: Edwards postulate 

The above Vogel-Pulcher law can also be understood using Edwards' analogy with thermodynamics. Assume again that all 
blocked states are equiprobable. Then the states that are most likely to be observed are the most numerous ones: this is the 
essence of statistical thermodynamics. Imagine for example that a container is divided by a movable piston, with grains in 
each compartments. The volume of the two compartments are Vi and V2, with Vi + V2 = V . The Edwards entropy <S is the 
logarithm of the number of blocked states for a given overall volume. The total entropy of the container is: 

5t =51+52. (119) 

The most probable value of V\ is such that this entropy is maximum, under the constraint Vi +V2 = V. Therefore, one has: 

In thermodynamics, this quantity (which is constant throughout the system) is equal to P/T, where P is the pressure. In the 
context of granular materials, dS/dV > was called the inverse compactivity by Edwards, and noted 1/X. 

A reasonable requirement for S{V) is that it should be proportional to the number of grains A'', and only depends on the 
free-volume fraction <i>: SiV) — Ns{^). One also expects that s{^) vanishes for "l? ^ 0. A possibility, suggested in [65,91], is 
that each grain has a number of possible positions proportional to the free-volume, which leads to s("i>) = log3>/vo, where vo 
is again a 'quantum' of free- volume. Using this form for s($), we finally find: 

showing that the Edwards 'temperature' vanishes for $ 07 Now, consider a small region immersed in a large container 
with grains, which acts as a reservoir of free volume. Exactly as in thermodynamics, one can show that the probability that a 
free-volume equal to v is found in this region is given by: 

Kv)ocexp(-J) (122) 

The probability for a hole of the size of a grain v to appear is therefore: p(v = w) ~ exp(— $*2/$), where $* = p*v. Assuming 
that the grain motion takes place when a suflBciently large hole appears, one finds that the rate of compaction has the same 
exponential form as above. 



Note the following intriguing consequence of Edwards postulate. Suppose that the two compartments contain grains of the 
very same material, but with different sizes (volume) ui and V2, such that Nivi — N2V2, i.e. the volume fraction occupied by 
the grains in the two compartments is the same. The equality Xi = X2 then imposes that $i/$2 = V2/V1, i-e. that smaller 
grains will have a higher free volume. This conclusion holds for more general choices of the entropy s{^). 
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C. Granular dynamics and the trap model 



The above 'free-volume' ideas are interesting and certainly contain important physical ideas. However, the description of 

the dynamics using a single macroscopic degree of freedom (namely the average free volume density) is too naive to account 
for the more sophisticated 'memory' experiments. Indeed, any rate equation of the form: 

^ = -7($,r)$, (123) 

is unable to explain why a packing prepared under different tapping conditions (different F), but such as to reach the same 
value $0 would evolve differently if tapped with the same amplitude. The same argument was used by Kovacs in the context of 
glasses to suggest that additional parameters are needed to fully describe a glassy state. Take the example of a one dimensional 
ferromagnetic Ising model that relaxes towards its equilibrium state a temperature T. Since the system does not order, the 
equilibrium state is characterized by a density of domain walls (or kinks), defining a characteristic domain size, which is the 
distance between the kinks. The non equilibrium state can also be characterized by a density of kinks, which decays with time 
towards the equilibrium. However, the full description of the non equilibrium state requires the specification of the distribution 
of the domain size - only its first moment is fixed by the average density of kinks. Different initial distributions of domain sizes 
with the same average value do evolve difi'erently at a given temperature. For example, if there is an initial excess probability 
of large domains (but such that the density of kinks is the equilibrium one), these will immediately break down, leading to a 
temporary increase of the density of kinks. 

A simple model that encapsulates both the spatial heterogeneity and the intermittency of the dynamics, is the 'trap' model. 
One should think of a glassy system as made of independent subsystems of a certain size ^ (see Fig. 10). Inside each of these 
regions, the dynamics is 'coherent' in the sense that hopping between different metastable states involves all particles within 
a blob of size < 5. Within each of these subunits, the dynamics can be thought of as a random walk in a rugged landscape, 
an idea with a long history in the context of glasses [74], and witnessing a strong recent revival in different contexts [85,37,75]. 
However, an element which often seems to be missing from the discussion is the fact that the dynamics of system as a whole 
necessarily results (for short range interactions) from the evolution in parallel of many subsystems: the dynamics of a particle 
in a given subregion is completely unaffected by the dynamics of fax away particles. An open problem is to identify precisely 
the (possibly time dependent) coherence length ^, and understand its temperature and density dependence (see [55] for very 
interesting results on this aspect in the context of Lennard- Jones systems) . 



FIG. 10. Schematic view of the dynamics in glassy material. Within each 'blob', the dynamics is 'coherent' in the sense 
that hopping between different metastable states involves all particles within a blob of size < ^, whereas the motion of far 
away particles does not affect the dynamics within this blob. Within each of these subunits, the dynamics can be thought of 
as a random walk in a rugged landscape, with an instantaneous (random) trapping time Tj in the ith blob. Grey regions are 
particularly 'slow' (r » tw)- 

The dynamics of a given subsystem can be seen as a succession of hops between different metastable states, each of which 
blocking the dynamics for a time that depends on the local packing fraction: high densities leads to long trapping times. The 
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state of the system at time t is described by a probability distribution P(4>,t), which counts the number of subsystems with a 
local value of the free volume. Each time a subunit unjams, we assume that it falls back into any of the blocked states with 
an a priori distribution Po(0), that reflects the number of states with a given packing density. The time evolution of P{(j),t) is 
then given by: 



This equation was introduced in [58,20,100] in the context of glasses, and extended to granular media by D. Head [77], where 
the free volume density </> plays the role of the energy E. This equation neglects any coupling between nearby sub-systems. 
This is obviously unrealistic since the free volume liberated at one point will in general help nearby volumes to unjam. A 
generalization of the above model that takes into account, to some extent, this coupling, can be found in [20,100,76]. 

The trap model, and its 'SGR' (Soft Glassy Materials) generalisation covered in Mike Gates lectures [134,66], exhibits a 
number of interesting features, such cis a genuine glass transition, aging and non linear rheology. A crucial ingredient of the 
trap model is the possibility for the average trapping time to diverge. Talcing t{4>) ~ ex-p{A/<p), where ^4 is a certain constant 
that may depend on the tapping amplitude F, and Po{4>) ~ cxp[^^s(</))] with s{<t>) = Blog(t>, we find that the distribution of 
trapping times "I'(r) decays very slowly, as l/r(logr)", where a = 2 + B£^'\ This means that one is always in the glassy phase 
of the trap model, since (t) = oo.^ Gorrespondingly, one expects not only slow (logarithmic) compaction but also aging effects, 
such as reported in numerical simulations [107]. In the glassy phase, the trap model describes the dynamics of each subsystem 
as essentially intermittent: either the system is blocked, or it moves fast to a quite different configuration; most of the time 
is spent in one particularly well jammed configuration. This intermittent dynamics of glassy systems begins to have some 
numerical [54,53,10] and experimental [120,32,39] support. However, one should remember that if one looks at the system as a 
whole, the activity will be dominated by the fastest regions (which can in turn, through the coupling between nearby that we 
have neglected above, unlock the slow regions). 

Interestingly, the trap model is able to reproduce some of the cycle effects reported above: the irreversibility in a tapping 
cycle [77], or the Kovacs effect [12,143]. Variants of the trap model can also be studied, where each 'trap' is decorated by 
the dynamics of smaller length scales [17,23] in order to account for memory and rejuvenation effects observed in spin-glasses 
and other glassy systems. The Sinai model is in this family of 'multi-scale' trap models [87,90,121], and was introduced in the 
context of granular media precisely to understand the memory and cycle effects [95], although the simpler 'monoscale' trap 
model seems to be able to account for most of them. It would be very interesting to exhibit experimentally a truly 'multiscale' 
dynamical phenomenon in granular materials. This would have the obvious advantage over many other glassy systems that the 
underlying mechanism, in terms of embedded length scales, can be directly observable [23,25]. 

There would be much more to say about glassy dynamics in the context of granular materials, but my lectures stopped at 
this point, after only quite simple ideas were expressed. Since some of the more elaborated concepts are common to spin-glass 
dynamics and glassy rheology, it is appropriate to refer to the lectures of Mike Gates, Leticia Cugliandolo and Giorgio Parisi 
in this volume, and also to [11,92,22,6] for further developments. 
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